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Abstract 


Variations of instantaneous heart rate appears regularly oscillatory in deeper levels of 
anesthesia and less regular in lighter levels of anesthesia. It is impossible to observe this 
“rhythmic-to-non-rhythmic” phenomenon from raw electrocardiography waveform in cur¬ 
rent standard anesthesia monitors. To explore the possible clinical value, I proposed the 
adaptive harmonic model, which fits the descriptive property in physiology, and provides 
adequate mathematical conditions for the quantification. Based on the adaptive har¬ 
monic model, multitaper Synchrosqueezing transform was used to provide time-varying 
power spectrum, which facilitates to compute the quantitative index; “Non-rhythmic- 
to-Rhythmic Ratio” index (NRR index). I then used a clinical database to analyze the 
behavior of NRR index and compare it with other standard indices of anesthetic depth. 
The positive statistical results suggest that NRR index provides addition clinical infor¬ 
mation regarding motor reaction, which aligns with current standard tools. Furthermore, 
the ability to indicates the noxious stimulation is an additional finding. Lastly, I have 
proposed an real-time interpolation scheme to contribute my study further as a clinical 
application. 


Keywords: instantaneous heart rate; rhythmic-to-non-rhythmic; Synchrosqueezing trans¬ 
form; time-frequency analysis; time-varying power spectrum; depth of anesthesia; electro¬ 
cardiography 
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Chapter 1 


Introduction 


1.1 Field of Project 

The human body undergoes unique conditions when under anesthesia. Anesthesiologists 
administer anesthetics, also known as anesthetic medicine, dependent on the patient’s 
physiological information. The physiological information may include the patient’s heart 
rate, blood pressure, and electroencephalography (EEG) as this information reflects the 
inner dynamics of the human body. My study will present a new perspective on the 
physiological information that could contribute to the clinical practice of anesthesia. 


1.2 An Obscure Phenomenon 

As a practicing anesthesiologist, I have noticed that the instantaneous heart rate (IHR) 
fluctuates dependent on the various levels of anesthesia - it became regularly oscillatory in 
deeper level of anesthesia, whereas it became less regularly in lighter level, (see Eig.1.1) This 
phenomenon, I refer to as ^^rhythmic-to-non-rhythmic phenomenon'', has been described 
qualitatively in the past. However, to the best of my knowledge, no quantitative method 
has been proposed. Also, the underlying physiological mechanism has not been addressed. 


1.3 Quantification for Depth of Anesthesia 

This phenomenon has induced three questions - 1) how to quantify this phenomenon into 
a useful index, 2) whether the quantitative index can reflect the level of anesthesia, and 
3) what are the underlying physiological mechanisms. 
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Rhythmic 


Nonrhythmic 



ECG 







Figure 1.1; 

The tracing of R-R peak interval (RRI) represents instantaneous heart rate (IHR) ex¬ 
hibiting “rhythmic” or “non-rhythmic” patterns during anesthesia. This phenomenon is 
difficult to be read from the raw electrocardiography (ECG) waveform. Interestingly, the 
ECG waveform is mandatory in standard anesthesia monitor, but IHR is not. 


For an anesthesiologist, information from this phenomenon is unavailable in the clini¬ 
cal perspective. From Fig.1.1, it is apparent that the “rhythmic-to-non-rhythmic” phe¬ 
nomenon can be seen from the tracing of instantaneous heart rate, but not from raw 
electrocardiographic (ECG) waveform. Although real-time monitoring of ECG is manda¬ 
tory in the modern anesthesia monitor, the tracing of instantaneous heart rate is not dis¬ 
played. Even though an anesthesiologist was aware of this phenomenon, he cannot obtain 
“rhythmic-to-non-rhythmic” phenomenon information from modern anesthesia monitor to 
improve the anesthetic management. 

If NRR index can provide additional clinical value, I believe that it could help in future 
developments of standard anesthesia monitor. 


1.4 Project Motivation and Goals 

My clinical observation leads to an intuition that this phenomenon arises from the in¬ 
ner physiologic dynamics under surgery and anesthesia. Hence, the quantification could 
provide valuable information to help anesthesiologists in administering anesthetics. 

Pursuing these questions involves model construction based on knowledge of three distinct 
fields: signal processing, clinical anesthesiology and physiology. Then, I will propose 
a quantitative index, referred to as “non-rhythmic-to-rhythmic ratio” (NRR) to quantify 
the phenomenon. The computation of NRR index requires multitaper Synchrosqueezing 
transform, a recently developed time-frequency analysis technique. 
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I used a clinical database to explore the property of NRR index. Meanwhile, NRR index 
was compared with standard depth-of-anesthesia indices used in current clinical practicing 
to understand its performance. 


1.5 Thesis Outline 

My thesis will begin with a literature review of the clinical anesthesiology, neural physiol¬ 
ogy, and time-frequency analysis. Bringing multidisciplinary knowledge up to date leads to 
a theory for the possible mechanism of rhythmic-to-non-rhythmic phenomenon. I proposed 
a model as the foundation of further analysis using time-frequency analysis. The adaptive 
harmonic model is compatible in both ways; fitting the descriptive property in physiol¬ 
ogy, and providing adequate mathematical conditions for the quantification. Based on the 
adaptive harmonic model and multitaper Synchrosqueezing transform, NRR index was 
developed to quantify rhythmic-to-non-rhythmic phenomenon. I used a clinical database 
to analyze the behavior of NRR index, and compare it with other current standard indices 
of anesthetic depth. The statistical result shows positive value in clinical; NRR index 
provides information regarding motor reaction earlier and better than all the other depth- 
of-anesthesia indices. The result suggests a possible anatomic area and functional part 
in the brain that mediate the rhythmic-to-non-rhythmic phenomenon. To benefit each 
individual patient undergoing surgery and anesthesia, NRR index should be implemented 
in real-time. An optimal spline interpolation technique is proposed for the development 
of “real-time NRR index”. 



Chapter 2 


Background and Previous Work 


2.1 Historical Background in Anesthesia 

More than seventy years ago, Guedel made a remark on the connection between anesthesia 
and respiratory patterns[l]. It appears that there is a direct relation between anesthetic 
depth and respiratory patterns. A deeper level of anesthesia results in regular respiratory 
patterns while a lighter level of anesthesia results in irregular respiratory patterns. Dur¬ 
ing that era, depth of anesthesia were concluded from physical signs, such as the changes 
in pupil size and movement of the eyeballs. However, the introduction of muscle relax¬ 
ant in clinical anesthesia has hindered the observation of muscle movement and Guedel’s 
assessment method. 

Even though Guedel’s observations became less known, this “rhythmic-to-non-rhythmic” 
phenomenon has been documented in other literature. About twenty years ago, Kato used 
power spectrum to analyze IHR in anesthesia[2]. His study showed that the spectrum 
appears to be concentrated during deeper level of anesthesia; whereas the spectrum is 
dispersing in lighter level of anesthesia. Kato also mentioned the respiratory effect on 
this phenomenon. It is possible that the power spectrum technique helps to capture the 
obscure phenomenon in IHR at that era. 


2.2 Current View on Anesthesia 

Today, we would like to further explain this phenomenon through the advancement of 
various scientific fields in signal processing, anesthesiology, and neural physiology. 

Because of the advancement in anesthesiology, we now understand that anesthetics exert 
effect mainly on the brain. Distinct anatomic location mediates different functions, with 
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differential susceptibility to same anesthetics. Furthermore, each patient has his own 
unique genetic profile that renders diverse response to anesthetics. The inter-individual 
genetic variability, plus the dynamic influence of surgical procedure, make it difficult to 
determine appropriate anesthetic levels based on the dosage of anesthetics or single index 
of anesthetic depth. There is strong evidence that anesthetic management can improve the 
long-term health of patients, it is important to have monitoring instruments facilitating 
a comprehensive anesthetic management specific for each unique patient. Hence, the 
quantification of “rhythmic-to-nonrhythmic” phenomenon should be worthwhile. 

Current advancements in signal processing has resulted in the development of time-frequency 
analysis. In particular, reassignment and Synchrosqueezing technique were not readily ac¬ 
cessible in the past. Furthermore, current advancements in neural physiology has also 
provided more insight into the underlying mechanism of the “rhythmic-to-non-rhythmic” 
phenomenon, which has not been addressed according to my best knowledge. A more 
detail literature review regarding three main disciplines is as follows. 


2.3 Perspective from Anesthesiology 

Anesthesia comprises several components, including hypnosis, analgesia, immobility, am¬ 
nesia, and autonomic nerve system stability[3-6]. For each patient, the dosage for different 
anesthetics is tailored to the patient’s specific surgical needs. Therefore, anesthesiologists 
need to continuously monitor the patient’s response to the anesthesia. Bispectral Index® 
(BIS), an index computed from electroencephalography(EEG), is the current standard 
tools used to monitor depth-of-anesthesia for either the hypnotic component or the level 
of sleepiness. However, analgesia component, commonly known as the inhibition of pain 
sensation in anesthesia, requires an indirect interpretation of physiologic information such 
as heart rate (HR) and blood pressure [7-11]. It is impossible to achieve optimal anesthe¬ 
sia using only one kind of anesthetics. Also, anesthesiologists are unable to monitor the 
depth of anesthesia by using merely one index, consequently inadequately measuring anes¬ 
thetic effects. For a more comprehensive administration of anesthesia, it may be helpful 
to integrate different types of information to come to a decision. 

It is important to understand what is wrong with simply relying on BIS monitor. BIS is 
a processed EEG, which measures the activity of the frontal cortical area only. On the 
other hand, amnesia is the suppression of memory, which is governed by medial temporal 
cortical area. Analgesia is the suppression of pain sensation, which is governed by various 
subcortical areas. Immobility is the inability of movement which is typically controlled 
by the thalamus and spinal cord[7, 12, 13]. Lastly, the entire autonomic nerve system is 
controlled by the brainstem[14-16]. Although BIS index is the current gold standard of 
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anesthetic depth monitor, BIS sensor is unable to measure the activities from all above 
listed anatomical areas except the frontal cortex. Literature in recent years indicates that 
anesthetics exert differential effects on different brain areas[3, 7, 17]. Thus, there is a 
considerable risk measuring frontal cortex to conclude the anesthetic effect on other brain 
areas. For example, researchers have reported that the BIS monitor cannot reduce the rate 
of intra-operative awareness, which requires adequate memory function suppression[18]. 


2.4 Anesthesia and Long-Term Mortality 

Anesthesia is more than merely falling asleep. Evidences support the connection between 
anesthesia and long-term mortality. Too deep of hypnotic index (BIS index) for too long 
in surgery is associated with worse one-year mortality[19]. It has also been reported that 
different anesthetics cause different immune responses, which affect tumor metastasis[20]. 

Evidence suggests incorporating “stress reduction” strategy into anesthetic management 
leads to longer life expectancy. First, infants who needed surgery for their congenital 
heart disease had lower long term mortality rate if they received adequate analgesics [21]. 
Second, more adult patients who underwent major surgeries survived one year later if they 
receive adequate medication to stabilize autonomic nerve system[22, 23]. Lastly, when 
epidural anesthesia, a regional anesthesia technique, is combined with general anesthesia 
for patients undergoing major surgeries, studies have showed that epidural anesthesia 
reduces stress and pain, leads to faster recovery, and overall better outcome [24]. There 
are more compelling studies that support the benefit of anesthesia management on long 
term outcome. These evidences suggest that a comprehensive optimization of anesthesia 
is worthwhile. 

In summary, anesthesia comprises several components. A comprehensive anesthesia man¬ 
agement that covers several components could benefit patient for many years. Liteurature 
survey also brings hope that quantitative analysis of “rhythmic-to-non-rhythmic” phe¬ 
nomenon could be useful to monitor the anesthetic effect on subcortical regions. 


2.5 Perspective from Physiology 

From a physiological perspective, the underlying mechanism of the “rhythmic-to-non- 
rhythmic” phenomenon could be partially explained by breathing mechanism and car¬ 
diopulmonary coupling[14]. It is known that the neural respiratory control comprises of 
two systems, I) the involuntary automatic control system and 2) the voluntary control 
system. The involuntary automatic control system is controlled by the respiratory center 
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in the brainstem. The voluntary control system is controlled by the forebrain[25]. The 
anatomies of these two systems are distinct from one another. The respiratory neural 
signals from these two systems compete with each other but are both integrated at the 
spinal cord to control the respiratory motor neuron. During spontaneous respiration, the 
involuntary automatic respiratory pattern generated in the preBotzinger complex [25-29], 
is rhythmic. On the contrary, the breathing pattern of the voluntary respiratory motor 
control is non-rhythmic, which involves larger cerebral areas, including cortical processing 
and thalamic integration[25, 30-34]. 

The suppression effect of anesthetics on the forebrain, including the cortex and the tha¬ 
lamus, is stronger than that on the brainstem[35]. Meanwhile, studies have shown that 
the preBotzinger complex, the involuntary respiratory control center, is less susceptible 
to the anesthetics[36, 37]. These relations suggest that under deeper level of anesthesia, 
the non-rhythmic respiratory is more suppressed than the rhythmic respiratory center. 
Hence, the IHR exhibits oscillation that is more rhythmic during deeper anesthetic level. 
Contrarily, the non-rhythmic respiratory center, being less suppressed during lighter level 
of anesthesia, causes activity more non-rhythmic in IHR. This concludes that IHR exhibits 
more nonr-hythmic. In summary, literature in physiology supports the correlation that 
the differential anesthetic effects between the involuntary control system and voluntary 
control system induces the “rhythmic-to-non-rhythmic” phenomenon in IHR. 

Thus, I hypothesize that IHR exhibits the central respiratory activity via cardiopulmonary 
effect [14], and it reflects the integration of rhythmic and non-rhythmic respiratory activ¬ 
ities. Although more evidence is necessary to clarify this hypothesis, I propose using the 
NRR quantification methodology as a potential tool to evaluate the depth of anesthesia 
from a different aspect than using EEG-based monitoring. 


2.6 Physiology of Amplitude and Frequency Modulation 

The human body constantly regulates the respiration to meet the requirement of metabolism. 
In other words, the human body has to inhale oxygen (O 2 ), and expel carbon dioxide (CO 2 ) 
adequately. The respiratory control center in the brainstem receives the feedback infor¬ 
mation of oxygen concentration and the acidity to modulates the volume and the rate of 
breath. The acidity is due to the fact that CO 2 exists as hydrogen ion (H"*") and carbonic 
acid in human body[16, 38]. 

Since CO 2 are the metabolic products of human body, the H+ concentration and O 2 
concentration in the human body are gradually, but consistently, changing. Thus, we 
can describe that the respiration signal in the brainstem is oscillatory with amplitude 
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modulation (AM) and frequency modulation (FM). The AM and FM are continuous and 
slow change in time. These physiologic features can be depicted in mathematical language 
for subsequent modeling and quantification in the next chapter. 

In summary, literature from physiology supports an integrated rhythmic and non-rhythmic 
neural activity reflecting the level of anesthesia. From a clinical anesthesiology perspective, 
technological advancements on monitoring various aspects of anesthetic depth can help 
improve the quality of overall patient care. 

2.7 Instantaneous Heart Rate and Heart Rate Variability 



Figure 2.1: A schematic diagram for the R-R interval (RRI) signal derived from the 
raw lead III electrocardiographic (EGG) signal. We refer to RRI signal as instantaneous 
heart rate. The diagram also showed the EGG-derived respiration (EDR) 

The term heart rate is important to all anesthesiologists and is usually referred to as the 
number of heartbeats within one minute. To an anesthesiologist, an ongoing heart rate 
indicates the stability of the heart condition, or the stress responses to surgical stimulation. 
The heart rate is mainly controlled by both sympathetic nerve and parasympathetic nerve 
system from the brainstem to the heart[16]. Due to the dynamic nature of patient’s under 
surgery and anesthesia, we do not feel comfortable with the heart rate information at 
an one-minute interval. Instead, a standard requirement in anesthesia practice includes 
a continuous real-time display of the heart rate. In this thesis, the “rhythmic-to-non- 
rhythmic” phenomenon is only observable at the instantaneous level. Conversely, the 
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average heart rate of a normal routine excludes the information of “rhythmic-to-non- 
rhythmic” phenomenon. 

Both the Instantaneous heart rate (IHR) and ECG derived respiration (EDR),(Fig.2.1) can 
be derived from the raw ECG waveform[39-41]. Although we can say that the “rhythmic- 
to-non-rhythmic” phenomenon in IHR is a kind of “generalized” heart rate variability 
(HRV), a medical guideline has set up the standard method of HRV analysis[42]. Un¬ 
fortunately, the standard analysis techniques stated in this guideline cannot reveal the 
information of the “rhythmic-to-non-rhythmic” phenomenon. In particular, the SDNNx, 
pNNx or the power spectrum technique measures information in a specified period, which 
inevitably neglects the dynamic and transient changes within this period. The guideline 
discusses but does not resolve the potential issue of stationarity in the data. This also 
supports my desire to rely on the current scientific advancements on signal processes to 
resolve the current medical shortcomings in my project. 


2.8 Integration of Multidisciplinary Backgrounds 

Integration of multidisciplinary backgrounds provides partially theoretical background for 
the mechanism of “rhythmic-to-non-rhythmic” phenomenon. Below is a list of key facts 
from various disciplines; 

1. In respiratory control, non-rhythmic activity originates from the forebrain, whereas 
rhythmic activity originates from the brainstem. 

2. The forebrain has far more neurons and synapses when compared to the brainstem. 

3. Neural structures with more synapses are more susceptible to anesthetics 

4. Respiration have an coupling effect on IHR 

Fact I and 2 are knowledge in neurology and neural anatomy. Fact 3 is from anesthesiology, 
and fact 4 is physiology knowledge. Each single fact seems ordinary when considered in 
its own disciplines. However, when I combine them together, the theoretical ground is 
justified. Then, we can proceed to the modeling problem in the next chapter. 



Chapter 3 


Theory and Modeling 


3.1 Proposal 

The “rhythmic-to-non-rhythmic” phenomenon in IHR[2, 43] suggests that IHR during 
anesthesia consists of two components: one more rhythmic, less affected by anesthesia, 
and one more non-rhythmic, suppressed by deep anesthesia, and the strength of these two 
components vary according to the anesthesia depth. 

I hypothesize that quantifying the ” rhythmic-to-non-rhythmic” phenomenon will reflect 
the level of anesthesia. In order to quantify this change in IHR, I will introduce a novel 
index which I referred to as Non-rhythmic to Rhythmic Ratio (NRR) index. 


3.2 Description of Proposed Methodology 

My hypothesis will be tested through a four step method. I will capture the changes in 
rhythm through IHR using the adaptive harmonic model. Next, I will use the multitaper 
Synchrosqueezing transform to extract the dynamic behavior from the ECG signals. In 
order to quantify “rhythmic-to-non-rhythmic” phenomenon, I will use the NRR index. 
Lastly, I will use a clinical database to compare the NRR index with other anesthetic 
depth indices. 


3.3 Adaptive Harmonic Model 


It is important to consider which appropriate model to utilize to adequately describe 
and quantify this rhythmic-to-non-rhythmic phenomenon. Specifically, it is important 


10 



Chapter 3. Modeling 


11 


to recall the high difficulty on describing the inner dynamic of the human body through 
mathematical formulas. This difficulty is due to the fact that the human body is a complex 
system creating non-linear and unpredictable interactions with the outer environment. The 
drastic influences of surgery and anesthetic medication complicates the difficulty further. 
It is fairly simple to ruin the capability of a model to describe the reality by imposing 
unrealistic conditions. That is the most challenging task currently at hand when we try to 
model the biological phenomenon of our study. In order to do so we have to be extremely 
cautious on what conditions we choose and impose on our study model. 

After careful observation from a time-varying spectrum, I propose a relatively conservative 
model to analyze the rhythmic part of IHR. As is shown in Fig. 1.1, the definition of 
rhythmic implies that the oscillation exhibits time-varying frequency and time-varying 
amplitude. However, the modulations are in a relatively slow rate. This presents the notion 
that the relative predictable oscillatory activity can be perceived as a rhythmic movement. 
Hence, a more conservative phenomenological modeling, also known as adaptive harmonic 
model, will be proposed as follows. 

The adaptive model integrates with the Synchrosqueezing transform creating a whole new 
scheme starting from the theoretical grounds to the clinical phenomenon. We start from 
motivating the technical model we consider to analyze IHR(t). It is commonly observed in 
clinical settings that stronger anesthetic levels are associated with more regular patterns 
in the instantaneous heart rate IHR(t), which is caused by the cardiopulmonary coupling 
effect under anesthesia [2, 43]. Based on this major observation, our treatment of the 
IHR(t) signal will be purely phenomenological; that is, the parameters and indices we will 
derive from the IHR(t) signal will be based solely on these signals themselves, and not 
on explicit, quantitative models of the underlying mechanisms. Since we mainly extract 
IHR from ECG signal as RRI, we will consider the following phenomenological model for 
RRI(t): 

RRI(t) = T{t) + A{t) cos(27r(/)(t)), (3.1) 

where T{t) > 0, A{t) > 0 and > 0. We call T{t) the trend, which captures the 
“average heart rate” over a long period; we require the trend to be positive, but not to 
be constant; we allow it to vary in time, as long as the variations are small enough. We 
call the derivative (f)'{t) of the function (l){t) the instantaneous frequency (IF) of RRI(t); 
we require IF to be positive, but not to be constant; we allow it to vary in time, as long 
as the variations are slight from one period to the next, i.e. |i;/>^^(t)| < ecj)'{t) for all time t, 
where e is some small, pre-assigned positive number. Likewise, We call A{t) the amplitude 
modulation (AM) of R{t), which should be positive and smaller than T{t), but is allowed to 
vary slightly as well, i.e. |A'(t)| < ecj)' {t) for all time t. In summary, we have the following 
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conditions for the model 3.1 [44, 45]; 

0 < A{t) < T(t), > 0, (3.2) 

f all t. (3.3) 

If a RRI signal satisfied the model (3.1), we call it rhythmic] otherwise we call it non¬ 
rhythmic. It is possible that a function can be composed of both the rhythmic and non¬ 
rhythmic components. 

This seeming complicated model is actually a direct generalization of the Fourier harmonic 
model [44, 45]. Close inspecting the waveform of IHR, we can see that the conditions on 
(l)"{t) and l^^(t)l are reasonable. In other words, the rhythmic RRI exhibits an oscil¬ 
lation with “small enough” frequency modulation and amplitude modulation. Without 
adequate knowledge of underlying mechanism, the adaptive harmonic model fulfills the 
phenomenological behavior of rhythmicity feature. 


3.4 Non-rhythmic to Rhythmic Ratio 

The amount of variance inside the RRI signal are known as Heart Rate Variability (HRV). 
HRV has been proven to be closely related to the autonomic activity [42]. The most 
commonly used tool to quantify HRV is the power spectrum (PS) [42]. As useful as the 
PS is, however, in this study the PS cannot capture the momentary dynamics in the RRI 
time series. This issue limits the application to human under highly dynamic situation 
~ in particular, the “rhythmic-to-non-rhythmic” phenomenon in anesthesia. The main 
technical focus of this paper is to resolve this limitation. 

In this paper I will quantify the rhythmic-to-non-rhythmic pattern of the RRI series by us¬ 
ing a novel index referred to as Non-rhythmic to Rhythmic Ratio (NRR). NRR is motivated 
by the clinical finding that deeper anesthetic levels are associated with the appearance of 
more regular oscillatory patterns in IHR. Thus, under deeper anesthesia, IHR oscillates 
more regularly. We call this regular oscillatory IHR “rhythmic”, which appears as a sharp 
peak on the PS. On the other hand, when the subject is under lighter anesthesia, IHR 
varies irregularly, and we call this IHR “non-rhythmic”. The non-rhythmic IHR exhibits 
irregular and random-like behavior that appears as a plateau on the PS. As a result from 
this observation, the RRI time series is composed of a non-rhythmic component and a 
rhythmic component with various ratios. Based on the above facts, we hypothesize that 
the time-varying ratio change of the non-rhythmic and rhythmic components may reflect 
the anesthetic effect on the brain. 
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To quantify the ratio of the rhythmic component and non-rhythmic component, we will use 
a new signal processing tool to fully capture the momentary RRI time series change. We 
will replace the PS with the time-varying power spectrum (tvPS) by applying a recently 
developed time-frequency analysis technique [43, 46], which is referred to as multitaper 
Synchrosqueezing transform [44, 45, 47]. The tvPS is a non-parametric generalization 
of the PS providing instantaneous PS information. We will then be able to trace the 
momentary strength of the rhythmic component and non-rhythmic component through 
analyzing the tvPS. 

By using the tvPS we will further extend the definition of the classical frequency do¬ 
mains HRV parameters. Recall that traditionally, the high frequency (HP) power, the low 
frequency (LF) power, and the low frequency to high frequency power ratio (LHR), are 
determined from the PS of the RRI time series. Vagal activity mainly contributes to HF 
power, meanwhile, sympathetic activity influences both LF power and LHR [42, 48]. By 
using the tvPS of RRI, we are able to obtain time-varying HRV parameters, including the 
time-varying high frequency power (tvHF), the time-varying low frequency power (tvLF) 
and the time-varying low frequency to high frequency power ratio (tvLHR). All these pa¬ 
rameters are suitable to provide a better understanding of the dynamical situations of 
anesthesia. 

The main index NRR is defined as the ratio between momentary non-rhythmic power to 
momentary rhythmic power. In mathematical terms, we compute the rhythmic component 
power by identifying the location of peak on tvPS. The non-rhythmic power is then defined 
as tvHF subtracted by the momentary rhythmic power. Finally, NRR is defined as 

, /Non-rhythmic power\ 

NRR = logio .. ■ - > 

\ Rhythmic power / 

where the value of NRR is high with non-rhythmic RRI and the value of NRR is low with 
rhythmic RRI. As previously discussed, the dynamical behavior of RRI(t) is captured by 
tvPS. Please see Figure 3.1 

In practical calculation, the non-rhythmic power is calculated as the rhythmic power sub¬ 
tracted from the high frequency power of RRI. The rhythmic component power is calcu¬ 
lated as the sum of the power near the IF of the rhythmic component within the range 
of high frequency on the tvPS. We define the IF of the rhythmic component, denoted as 
fr{t), as the time varying position of the maximal power within the high frequency range 
on the tvPS of RRI(f): 


fr{t) := argmaxMRRpx(L6> 


(3.4) 
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Figure 3.1: Fifty-second electrocardiography (ECG) waveforms, time tracings of R-R 
interval (RRI), ECG-derived respiration (EDR), and the corresponding spectra in lighter 
level of anesthesia (A, B, G) and deeper level of anesthesia (D, E, F). Whether in lighter 
(A) or deeper level of anesthesia (D), it is hard to read the IHR information from raw 
EGG waveform by naked eyes. Time tracings (B, E) of RRI and EDR and their power 
spectra (G, F) show that RRI is more non-rhythmic in lighter level of anesthesia then in 
deeper level. The frequency bands around the instantaneous frequency /r and its multiple 
frequencies (the unshaded area) constitute the rhythmic component. The rest (shaded 
area) constitutes the non-rhythmic component. dB = decibel 
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which is implemented by the curve extraction algorithm considered in [49]. Next, we will 
calculate the rhythmic component power, denoted as Pr{t), as the sum of the power inside 
the bands around fr{t) on the tvPS of RRI(t). In particular, 

l'nfr(t)+0.01Hz 

Prit) ■■= / 

J nfr{t)—0.01Hz 

Here the width of the band is chosen to be O.OlHz in an ad hoc way. We mention that when 
the signal is of broadband spectrum, there is no dominant curve that can be visualized, 
as shown in Fig.3.1. In this case, the most dominant curve defined in (3.4) from the 
time-frequency representation is viewed as the rhythmic component. 

The non-rhythmic component power is defined as the power within the high frequency 
subtracted by the rhythmic component power. Lastly, NRR is defined intuitively as the 
ratio of the non-rhythmic component power to the rhythmic component power: 

Through defining power within absolute boundary as Pr{t), NRR can differentiate the level 
between “rhythmic” and “non-rhythmic”. In this thesis, the variables affixing (t) mean 
that they are functions of time. It is important to notice that these sorts of time-varying 
parameters describe the momentary behavior of the signal. Since this instantaneous in¬ 
formation cannot be captured by the classical power spectral analysis, we would expect to 
gain more insight into the IHR signal by considering our definitions. 


3.5 Non-rhythmic Component vs. Stochastic Process 


For the purpose of this thesis, the term non-rhythmic is defined as any signal that appears 
“irregular”, “noisy”, or “as random”. This poses two interesting questions. First, does 
this mean that we should treat all non-rhythmic signals as a stochastic process? Second, 
will we be able to model the non-rhythmic signal as an auto-regressive moving-average 
(ARMA) model or other stochastic process in order to analyze it using parametric meth¬ 
ods? However, I would like to also remind that we should consider these issues using our 
medical knowledge. 

One of my hypotheses in the present study is that the non-rhythmic signal reflects highly 
structured and complex activity of the human brain. Although it is too early to determine 
that the human brain is a “deterministic” or “stochastic” system, when awake, the human 
brain can retain memory for longer periods of time to execute complex tasks. The thought 
or idea in the brain lasts longer too. Therefore, the brain exhibits non-rhythmic activity 
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when its inner state is longer lasting. However, the above statements would contradict the 
stochastic process behavior that presents no memory effect in a Gaussian random process. 
That is, the autocorrelation function is a Dirac’s delta function: 


= ^yin)yin-l), 
i&z 

E[i?i?yy(0)] = a and K[RRyy{k)] = 0 for /c > 0 when y is a Gaussian process. 

Although non-rhythmic signals may appear similar to a stochastic process, the human 
brain that shows non-rhythmic activities should result in a more structured state. There¬ 
fore, a non-rhythmic signal in its natural state is not completely “stochastic”. This is why 
I would like to use the term ” non-rhythmic”. 

My second hypothesis is that there are two components, non-rhythmic and rhythmic, 
that are contained in the ” rhythmic-to-non-rhythmic” phenomenon. It is unknown how 
these two components coexist but one simple way to model their co-existence would be 
superposition. Regardless if we fully understand how the two components coexist, the fact 
that stochastic process is fundamentally a “white-noise-driven” process means it doesn’t 
match the “rhythmic component plus non-rhythmic component” concept. 

The above technical issues regarding parametric modeling for “rhythmic-to-non-rhythmic” 
phenomenon needs to be carefully addressed in the future work. To narrow down the 
problems I have to solve in my thesis, I decided to choose a non-parametric method as 
a starting point. Still, statistical signal processing technique is a viable direction for my 
future studies. 


3.6 Summary 

In summary, the proposed adaptive harmonic model (3.2) describes a fundamental physi¬ 
ologic property: the rhythmic component possesses frequency modulation and amplitude 
modulation. However, the law of physiology and chemistry sets a limit to these modula¬ 
tions. 

From the previously mentioned theory to the modeling and quantification of the NRR 
index, I keep seeking the time-varying ability to capture the momentary dynamics, which 
is essential to the dynamic nature of clinical anesthesia. This feature is also fundamentally 
distinct to theories of other signal analysis methods. After establishing the theoretic part, 
I will focus on the next problem: the time-varying spectrum. 



Chapter 4 


Time-Frequency Analysis, 
Reassignment, and 
Synchrosqueezing 

4.1 Time-Frequency Analysis 

Chapter 3 provided strong evidence that in order to execute the quantification of NRR 
index within anesthesia the adaptive harmonic model requires a strong time-varying spec¬ 
trum. A recently developed time-frequency analysis technique is called the multitaper 
Synchrosqueezing transform. I believe that this technique is the most useful for my pro¬ 
posed study. 

Time-frequency analysis is a type of signal processing tool that generalizes the power 
spectrum to provide time-varying ability. This type of analysis generates a series of power 
spectra that has the ability to represent instantaneous spectra. Short-time Fourier trans¬ 
form (STFT) and continuous wavelet transform (CWT) are two types of well known tools 
for biomedical signal analysis[43, 50-53]. Both of these tools are able to capture the instan¬ 
taneous frequency information of the biological signal. Huang and Chan were the first to 
employ time-frequency analysis to analyze the dynamic changes in IHR during anesthesia 
[46]. 

Both STFT and CWT require window functions to “focus” or “localize” the signal. Take 
STFT as example: 


/*oo 

vP{t,i) = / X{s)h{t - 

J —OO 


(4.1) 
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where /i is a window function chosen by the user, for example, the Gaussian function 
h{t) = Note that X{T)h{t — r) is the piece of signal around t chosen by the window 
function. We then define the spectrogram of X(t) [50] as 

= (4.2) 

From equation 4.1, it is apparent that choosing different window functions causes different 
results. The influence of the chosen window function in relation to the STFT is no less 
than the signal per se. In particular, the window function imposes a resolution limit which 
relates to the well-known Heisenberg’s uncertainty principal. The larger the time window 
h is, the higher the spectrum resolution is, but with low temporal resolution, and vise 
versa. 

The window function also influences other TF methods, such as CWT, and Smoothed 
Pseudo Wigner-Ville Distribution, to name a few. Rather than fiddling with various TF 
methods, various window functions, and different sizes, it is highly advised to choose a 
technique that alleviates the above issues. 


4.2 Reassignment 

To cope with the resolution limitation of time-frequency analysis, it is possible to further 
use the phase information which is discarded in the classical STFT spectrogram 

(4.2) and scalogram[50, 51]. Kodera et al. first proposed this technique[54] by calculating 
the instantaneous frequency and group delay Auger and Flandrin coined 

this technique as the reassignment and further generalized it to “Fourier transform based” 
time-frequency analysis[50]. The reassignment technique provides a sharper representa¬ 
tion, which improves the readability of the STFT spectrogram. Reassignment is helpful 
in visualizing the dynamics of the biological signal that contains multiple nonstationary 
components in the TF plane (Fig. 4.1). 

The phase information hidden in STFT contains information that explains why time- 
frequency domains are spreaded out. [55]. To further eliminate this issue, we can shift the 
STFT coefficients accordingly to the reassignment rules determined by the phase informa¬ 
tion. The temporal reassignment rule and frequency reassignment rule are as follows: 



Frequency (Hz) Frequency (Hz) 
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Figure 4.1: A side-by-side comparison between STFT spectrogram (A) and Reassign¬ 
ment spectrogram (B). The tracing of original R-R interval signal recorded in anesthesia 

is denoted as red line. 
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t + 3? 


vP\t,C) \ 

I 

vP{t,f) J 


where 3? denotes the real part and $5 denotes the imaginary part and h, Th, and Dh 
are the three related window functions. Here Th is the time-ramped window used to 
determine the temporal reassignment rule and Dh is the differential window used to de¬ 
termine the frequency reassignment rule. The reassignment is then achieved by applying 
the reassignment rules: 

/ CO POO 

/ -ts,^)6{r] - Z,i)^sdr], 

OO J —CO 

where we move the spectrogram coefficient V^\s,if} at time s and frequency ry to a new 
location and according to the reassignment rules. 

The reassigned spectrogram is a valuable tool for signals that contain dynamic information 
of anesthesia in my current study[43]. However, aside from the visual information repre¬ 
sented in the TF plane, it is also desirable to execute further analysis of the signals, such 
as isolating a component in TF plane, or performing reconstruction in time domain. Syn¬ 
chrosqueezing transform is a recently developed reassignment technique that can provide 
this possibility. 


4.3 Synchrosqueezing Transform 

Synchrosqueezing transform is a special kind of reassignment technique[44, 45, 56]. This 
technique uses a light-weight mathematical model (4.5). The light-weight mathematical 
model helps with further processing of signals, such as isolating a component in the TF 
plane or reconstruction in time domain. 

For my present study. Synchrosqueezing transform can be used to present non-parametric 
generalizations of the power spectrum. The Synchrosqueezing transform is able to provide 
the time-varying power spectrum (tvPS) which reveals the instantaneous changes that 
happen in IHR. In addition, comparing with traditional STFT and wavelet, it is visually 
more informative, more resistant to several kind of noise, and more resistant to the influ¬ 
ence of chosen window function [44, 45, 49]. Mathematical works have supported these 
abilities, as well as the ability of inverse transform. Luckily we do not need to impose too 
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Figure 4.2: A side-by-side comparison between STFT spectrogram (A) and Reassign¬ 
ment spectrogram (B). The tracing of original R-R interval signal recorded in anesthesia 

is superimposed as red line in panel B 


much effort on choosing parameters as using the Synchrosqueezing transform technique 
can easily capture the oscillatory activity. 

The detailed information on the algorithm leading to tvPS is as follows. It has been known 
that the phase information hidden in STFT contains more information that explains why 
the STFT spectrogram is spreaded out. A solution that could help remedy this issue is by 
shifting the STFT coefficients according to reassignment rules determined by the phase 
information. The frequency reassignment rule is as follows: 
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-idtVPihf) 


(4.3) 


where dt means the partial derivative with related to t. The Synchrosqueezing transform 
is then achieved by applying the reassignment rule: 

poo poo 

Sx\t,0-= / yx\t,r])5{\f-n{t,r])\)dr], 

J —oo J—oo 

where for each fixed time t, we move the STFT coefficient {t, rj) at frequency rj to a 
new location according to the reassignment rules 

One way that the Synchrosqueezing transform is unique to the reassignment technique 
is that the Synchrosqueezing transform has the ability to reconstruct back to the time 
domain. Suppose the signal X(t) analyzed by Synchrosqueezing transform is: 


X{t) = T{t) + A{t) cos{2tt (4-4) 
where T{t) > 0, A{t) > 0 and if'it) > 0- 

In other words, X{t) comprises of an oscillatory component and a trend component. 

In addition, the amplitude modulation and frequency modulation of the oscillatory com¬ 
ponent is under the condition as follows: 

0<A{t)<T{t),O{t)>0, (4.5) 

|^^(i)| < \0'{t)\ < ecj)'{t), \T{ty\ < e for all t. (4-6) 


Relevant results of the Synchrosqueezing transform can be found in literatures[49, 57- 
60]. A recent review article discussed in depth on the relation of the TF analysis, the 
reassignment and Synchrosqueezing technique [61]. 


4.4 Multitaper Estimation 

From a non-parametric spectral estimation perspective, Thompson proposed a powerful 
multi-window technique that improves spectral estimation of the random process [62, 63]. 
Thompson uses the discrete prolate Spheroidal sequence to define a series of window func¬ 
tion, which are orthogonal to each other, that have the same length in time domain. 
The average of the spectra from multiple windows reduces the spectrum variance but 
also prevents the bias of spectral estimation induced by varying window size. One strong 
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advantage of Thompson’s spectrum is that we are able to estimate the spectrum of a 
stationary random process without resort to parametric estimations. 

The principle of Thompson’s technique has been extended to the field of nonstationary 
signal. Multitaper in time-frequency analysis is designed to be used with Hermite functions 
as the sequence of multiple windows that are orthogonal to each other in time-frequency 
plane. The final spectrogram estimation can then be obtained from the average of multiple 
realizations. Taking the STFT spectrogram as an example, the multitaper spectrogram 
can be defined as follows; 


Qx,K{t,C) 


K 




K 


k=l 


The multitaper spectrogram and multitaper scalogram were shown to make success for 
nonstationary stochastic signal[64]. To achieve both high resolution for deterministic non¬ 
stationary signal and low variance for stochastic nonstationary signal, Xiao and Flandrin 
combined the reassignment technique and multitaper technique[47]. 

From the background information mentioned in Chapter 2 and the model in Chapter 3, it 
can be seen that the IHR under anesthesia should contains both non-stationary component 
and stochastic component. We have reported the advantage of using multitaper time- 
frequency reassignment in extracting these dynamic features [43]. 


4.5 Multitaper Synchrosqueezing Spectrogram 


Synchrosqueezing transform responds well to signals from different types of stochastic 
process. One dilemma we face is understanding the true mechanism in biological signal. 
Since the quantification of the non-rhythmic component is an important part of our current 
model(3.4), we have decided to apply the multitaper method to further optimize the 
performance of the Synchrosqueezing transform X{t). We have chosen the multitaper 
Synchrosqueezing spectrogram to achieve the tvPS as follows: 


Mx,K{t,C} 




k=l 


where K > 0 is the number of windows chosen. We choose the the k-th Hermite function 
as hk- For my current study, I will use RT = 10 and the first ten Hermite window functions 
will be hk, /c = 1,..., 10 when deriving tvPS. 
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Figure 4.3: The time-varying power spectrum of instantaneous heart rate (IHR) as data 
used in figure 4.2, realized by multitaper Synchrosqueezing spectrogram, showing rich 
dynamic features of human body in anesthesia. As shown in figure 4.2, the readability is 
better than STFT spectrogram, and also better than the reassigned spectrogram. 


The advantage of using a multitaper Synchrosqueezing spectrogram is summarized as 
follows: 


1. Robust to several different kinds of noise, which might be slightly nonstationary. 

2. Visually informative. When the signal is rhythmic, that is, satisfying the model (3.1), 
a dominant curve following the IF can be seen on the multitaper Synchrosqueezing 
spectrogram; otherwise the multitaper Synchrosqueezing spectrogram is blurred. 

3. Adaptive to the data so that its dependence on the chosen window function is weak. 

We compared various methods of TF analysis in Tabled. 1 and in simulated data analy- 
sis(Fig.4.4). 


4.6 Adaptive Filter and Parametric Estimation 

Another kind of techniques based on parametric modeling for time series data analysis 
is the adaptive filter, including the famous Kalman filter, and the parametric spectral 
estimation[65]. The adaptive filter is designed to use a cost function to represent the 
weighted mean square error. The goal is finding a set of parameters a* to minimize the 
cost function. 


Qi = argminE{e^[n]} 
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Figure 4.4: An artificial data, as a summation of two time-varying components and 
a Gaussian white noise, and its TF analysis realized by STFT spectrogram, Syn- 
chrosqueezed spectrogram and the multitaper Synchrosqueezed spectrogram. From 125 s 
to 250 s, it is apparent that Synchrosqueezing technique improves the readability of the 
two components, and multitaper technique refines further. 


If we have adequate knowledge about the signal, we can obtain the optimal point on the 
performance surface of the cost function by using partial differential. 

^E{e^[n]} = 0 
oai 

It is apparent that TF reassignment is fundamentally different from adaptive filter. TF 
reassignment is a non-parametric deterministic signal processing technique, whereas the 
adaptive filter is a parametric technique for stochastic signal processing. Regardless of the 
difference, both techniques use the same mathematical tools to achieve their goals. 

Kalman filter in essence is a real-time prediction-correction algorithm that provides para¬ 
metric estimation in the sense of weighted least-square fitting[62, 65]. As the physiologic 
mechanism underlying rhythmic-to-non-rhythmic phenomenon is not adequately under¬ 
stood, a technique which requires both elaborate model construction and careful parame¬ 
ter selection is not my first-line choice for my research problem. Besides, From Chapter 3, 
it shows that the NRR index has a high relation to oscillation. Therefore, for my study, 
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Table 4.1: Summary of abilities of TF methods. Details can be found in references cited 
in this chapter. X: unable; V: feasible; ?: questionable 



Reassignment 

SST 

STFT 

CWT 

Instantaneous frequency estimation 

good 

good 

fair 

fair 

Visual readability 

good 

good 

fair 

fair 

Component reconstruction 

X 

V 

V 

V 

Robust to noise 

? 

good 

fair 

fair 


the output of the Kalman filter as a matrix of state variables may not be as intuitive as 
the spectrum technique for data exploration and interpretation at this stage. 

Parametric spectral estimation has a complete theory of modeling based on stochastic 
process [62]. One critical issue of the parametric estimation is that the mechanism of the 
signal we are facing in the present study is unknown, and is constantly changing. The para¬ 
metric spectral estimation has a lower tolerance to incorrect modeling. In other words, it is 
crucial to ensure that the correct model are constructed to generate reasonable parametric 
spectral estimation. In addition, parametric spectral estimation does not behave well when 
the signal is a superposition of multiple components. The rigidity of model construction 
contrasts well with the flexibility of Synchrosqueezing: flexible in window selection and 
adaptive to signal. Certainly, it is possible that the technique realizing time-varying power 
spectrum in the future may belong to the family of parametric spectrum. 


4.7 Common Condition in Physiology and Mathematics 

The most unique characteristic of the Synchrosqueezing transform that is relevant to this 
project is the mathematical condition 3.2 which shows that the frequecy modulation and 
amplitude modulation are limited. The condition used on the Synchrosqueezing transform 
is exactly the same condition as the adaptive harmonic model in order to describe this 
rhythmic physiological behavior. Therefore, the Synchrosqueezing transform matches per¬ 
fectly with the adaptive harmonic model for us to quantify the “rhythmic-to-non-rhythmic” 
phenomenon that was previously discussed. 



Chapter 5 


Performance Evaluation using 
Clinical Data 

5.1 Clinical Database 

In my hopes to investigate the clinical value of NRR index, I have collected and analyzed 
a clinical database. The prospective observational study spanned the whole period of 
typical surgical anesthesia. ECG, BIS index, and the end-tidal concentration of sevoflurane 
were all recorded continuously. Anesthetic events, including loss of consciousness (LOG), 
skin incision, first reaction of motor movement during emergence period, and return of 
consciousness (ROG), were registered with precise time stamps. 

After receiving approval from the Institutional Review Board (Taipei Veterans General 
Hospital, Taipei, Taiwan), we enrolled 31 female patients. All patients were American 
Society of Anesthesiologists (ASA) physical status of I or II and were scheduled for la¬ 
paroscopic gynecological surgery undergoing general anesthesia. We obtained informed 
consent in paper from each of our patients. All ASA I or II status patients are generally 
under healthy conditions. It is important to investigate the changes of the NRR index 
under the normal physiological conditions. We are aware that some diseases can affect 
HRV (and IHR) untowardly. Thus, patients that are excluded were those younger than 20 
years old or older than 55 years old, have a body mass index below 19 or higher than 30 
and have either cardiovascular disease, arrhythmia, intake of medication that affects the 
neurological and cardiovascular system, potential cancer, diabetes mellitus, or anticipated 
difficult airway. All these surgeries were performed between the hours of 8 a.m. and 7 
p.m. 
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5.2 Anesthesia Protocol 

The standard anesthetic monitoring (HP agilent patient monitoring system) was applied to 
all patients. This included ECG, pulse oximeter (Sp02) and non-invasive blood pressure. A 
BIS electrode (BIS-XP sensor) and the ECG recorder (MyECG E3-80; Micro-Star Int’l Co., 
New Taipei City, Taiwan) were applied to the patients prior to anesthetic induction. After 
the initial step of cleaning the skin with rubbing alcohol, the BIS sensor is applied to the 
right side of the forehead. As a standard practice, all anesthetic managements, including 
the administration of medications, and airway management were under the discretion of 
the anesthesiologist. After intravenous infusion of lactated Ringer’s solution begins, we 
started anesthetic induction with fentanyl 2.5 - 3 pg/kg. Next, right after preoxygenation, 
hypnosis was induced with propofol 2 ~ 2.5 mg/kg. LOG was then assessed by no response 
to verbal command. Cisatracurium O.I5mg/kg was used to facilitate tracheal intubation. 
Subsequently, mechanical ventilation was started from volume control mode with oxygen- 
gas-sevoflurane mixture as low flow anesthesia. The respiratory rate and tidal volume of 
the ventilator was adjusted to maintain an end-tidal carbon oxide (ETCO 2 ) that ranges 
between 35 and 40 mmHg, and to keep the peak airway pressure lower than 25 mm Hg. 

The laparoscopic skin incision was made by a surgical blade. The adequacy of anesthetic 
depth (during the whole period of surgery) was determined by the BIS index (< 60) and 
based on the anesthesiologist’s judgment. If the anesthesiologist determined an inadequate 
level of analgesia, a bolus dose of fentanyl would be given to the patient. Toward the end 
of the surgery, the peritoneal gas was deflated and the patient’s body position recovered 
from Trendelenburg’s position to level-supine position. As the wound closure began, the 
anesthetic gas was on a consistent decrease. Muscle relaxation was reversed with a com¬ 
bination of neostigmine 0.05 mg/kg and glycopyrrolate 0.01 mg/kg during controlled ven¬ 
tilation. Once the patient regained spontaneous breathing the controlled ventilation was 
halted. Patients that exhibited inadequate spontaneous breathing (ETCO 2 >50 mmHg 
or Sp02 <95%) were assisted with manual positive ventilation, otherwise the patients 
continued to breath spontaneously until regaining consciousness without the interference 
of positive-pressure ventilation.ECG and BIS are the main electrophysiological signals we 
analyze for this project. When we observed the emergence period, physical contact to 
the patient was avoided and minimized if necessary to prevent the interference on the 
sensors of EGG and BIS monitor. Any unnecessary interference to patient’s respiration 
from the manual (bag) ventilation was avoided also. When the patient showed adequate 
consistency in their spontaneous breathing, we removed the endotracheal tube. If in any 
case there was inadequate ventilation or upper airway obstruction, corrective actions were 
taken that include either mask ventilation or nasopharyngeal airway insertion. During the 
emergence period, I carefully assessed first reaction as any first visible motor reaction such 
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as movement of the arms or legs, coughing, or grimace. ROC was assessed through the 
patient’s ability to open their eyes and follow simple procedural commands. 


5.3 Data Acquisition 

The BIS index was continuously recorded from an Aspect A-2000 BIS monitor (version 
XP, Host Rev:3.21, smoothing window: 15 seconds; Aspect Medical Systems, Nattick, 
CA, USA) connected to a laptop computer (Asus Corp., Taipei, Taiwan). During surgery, 
corrective measures were taken to improve signal quality of the BIS sensor when the signal 
quality index was lower than 50; otherwise the BIS data during this period was discarded. 
The raw limb lead ECG sampled at lOOOHz and 12-bit resolution, was recorded for off-line 
analysis. The clocks on the laptop and ECG recorder were synchronized with a time accu¬ 
racy of ±1 second. Time stamps of BIS record were provided by the laptop. The inhaled 
and end-tidal concentrations of anesthetic gas, which was sampled from the connection 
piece close to the endo-tracheal tube and detected by the gas analyzer on a Datex-Ohmeda 
S/5 anesthesia machine (GE Health Care, Helsinki, Finland), were also recorded on the 
laptop. All data was recorded without any interruption until the patient regained ROC. All 
data collected for analysis of sevoflurane concentration correlation and anesthetic events, 
including skin incision, first reaction and ROC, are under single ventilation mode. There 
were no transitions from the mechanical ventilation to the spontaneous ventilation or vice 
versa in the data. However, during LOG, the transition from the spontaneous ventilation 
to mechanical ventilation is possible. 

The offline ECG signal was analyzed in several steps. The R peak detection was auto¬ 
matically determined by taking lead I, H and HI ECG signals into account to improve 
accuracy. The ectopic beats were removed and interpolated from the data to eliminate 
incorrect R peaks and ectopic beats through visual verification. When electrocauteriza¬ 
tion is severely interfered by the ECG, the data segment were also discarded. Whenever 
electro-cauterization severely interfered the ECG, this data segment was discarded. IHR 
was then derived from the R peaks through the cubic spline interpolation. RRI were 
resampled to be equally spaced at 4 Hz from IHR for subsequent analyses (Fig.2.1). 

In order to correct the hysteresis between the end-tidal sevoflurane concentration and the 
anesthetic effect on the brain, the estimated effect-site sevoflurane concentration (Ceff) 
was derived from the end-tidal sevoflurane concentration (Cet) by the following first order 
differential equation. This equation is based on the pharmacokinetic-pharmacodynamic 
modeling [66]: 
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dCeff ^ 
dt 

where the constant K^q was assumed to 
/min according to previous studies [8]. 


Keo{Cet — Ces) 

be constant for all patients and defined as 0.20 


5.4 Statistical Analysis 

We treated BIS and ECG-derived continuous indices, including NRR, tvHF, tvLF, tvLHR 
and HR, as anesthetic depth indices. The performances of the indices were evaluated in two 
ways: the ability for these indices to predict anesthetic events and their correlations with 
effect-site sevoflurane concentration during emergence period. Because the concentration 
of anesthetic gas consistently decreased and the influence from surgery was relative minor, 
we considered the emergence period for correlation analysis. A p value less than 0.05 would 
be considered as statistically significant. Multiple significant tests were calculated by using 
Bonferroni correction. Statistical results are expressed as mean (standard deviation [SD]). 

Prediction probability (Pk) analysis is a versatile statistical method for measuring the 
performance of an anesthetic depth index [67]. A value of one means that the observation 
is always correctly predicted, a value of 0.5 indicates that the observation is predicted at 
a 50/50 chance. The ability to predict anesthetic events can be investigated by using the 
serial prediction probability analysis[67]. The correlation with sevoflurane concentration 
employed Pk analysis and Spearman rank correlation. Estimation of Pk and its standard 
error was obtained through the Jackknife method. Prior to using the null hypothesis, any 
Pk value less than 0.5 would be changed into one minus the Pk value[67]. 


5.5 Serial Prediction Probability Analysis in Predicting Anes¬ 
thetic Events 

To evaluate the performance of the anesthetic indices in predicting the anesthetic events, 
serial Pk analysis was designed as follows. The following timestamps are used as the 
baseline: one minute before LOG, five seconds before skin incision, three minutes before 
the first reaction and three minutes before ROG. The serial Pk analysis is performed 
as successive Pk analyzes of data pairs: index on baseline timestamp versus indices on 
subsequent successive timestamps spaced by five seconds. As a result, the out of the serial 
Pk analysis is a sequence of time-varying Pk values which reveals the temporal relation 
between anesthetic events and the indices. By plotting the serial Pk value and its standard 
error bar, it becomes possible for us to do hypothesis test simply with the naked eye. If 
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there is significant difference (p < 0.05), it can be established if 1.5 times of their standard 
error bars do not overlap [67, 68]. 


5.6 Correlations with Sevoflurane Concentration 

Correlation between the above indices and sevoflurane concentration was investigated in 
the emergence period. We chose the time interval where the sevoflurane concentration 
is monotonically and continuously decreased. The period of spontaneous breathing is 
defined from the start of adequate spontaneous breath to ROC. We employed the Pk 
analysis and Spearman rank correlation to analyze performance of indices sampled every 
4 seconds. The results were tabulated as weighted averages according to each patient’s 
data length. The bootstrap method was used to calculate the 95% confidence interval of 
Spearman correlation based on 10,000 samplings. A value of Spearman correlation closer 
to 1 is better since we expect that the indices will increase to correlate with the decrease 
of sevoflurane concentration. 


5.7 Algorithm of Serial Prediction Probability (sPk) Anal¬ 
ysis 

In this section we detail a variation of the Prediction probability (Pk) analysis [67], referred 
to as Serial Prediction Probability (sPk) analysis. sPk is employed in this study to evaluate 
the prediction performance of the BIS and ECG-derived continuous indices, including 
NRR(f), tvHF(t), tvLF(t), tvLHR(t) and HR. 

We consider K different measurement times, which are the times when we record the 
anesthetic indicators. Then we determine two special timestamps, one is referred to as 
the event time, indicating the happening of the event, for example, the first reaction, and 
one is referred to as the base time, which is T > 0 seconds before the event time. Here T 
depends on the event we are studying. Collect N subjects and record the measurement 
times and the event time of each individual. Note that the measurement times and the 
event time and are different among individuals, and we denote the k-th. measurement time 
of the n-th subject as k = 1,... ,K, and the base time as to,n- Denoted as the 
anesthetic indicator (e.g. BIS) of the n-th subject recorded at time tk^n, and denoted as 
the associated outcome (e.g. first reaction) at time tk,n- 

To evaluate the prediction power of an anesthetic indicator at the A:-th measurement time, 
we take the following two datasets. At the measurement time tk^n, we record the sam¬ 
ple anesthetic indicators {x^,..., x^} and the sample outcomes {Vi,... ,y%}, where each 
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pair is recorded from the n-th subject. Denote := {(xj^, yf),..., (x^, y^)} 

as the measurement data at the k-th. measurement time. Moreover, at the base time 
to,n, we record the sample anesthetic indicators {x®,..., x^^r} and the sample outcomes 
{y ^,..., y%}, where each pair is recorded at the base time. Denote Zq := {(xj*, yi), ■ ■ ■, (x^, 
as the data at the base time. 

For a given individual, there are five possible relationships between pairs of anesthetic 
indicator and outcome at the k-th measurement time and those at the base time; 

1 . (x^,y^) and {x^,y^) are concordant if they are rank ordered in the same direction; 

2 . (x^,y^) and {x^,y^) are disconcordant if they are rank ordered in the reverse di¬ 

rection. 

3- {Xn, Vn) and (x^, y^) tie in x if x^ = xj^ while y’f / yj^; 

4- {xi, yn) and (x^, y^) tie in y if y^ = y^ while x^ 7 ^ x^; 

5- {xi, yn) and (x^, y^) tie in both x and y if x^ = x^ and y^ = y^. 


Usually the anesthetic indicator is of finer scale and the outcome is of coarser scale. Thus, 
unlike the usual notion of correlation, in the Pk analysis the tie in x is undesirable but 
we tolerate the tie in y and the tie in both x and y. Based on this notion, we define 
the serial Pk index in the following way. Denote Pdk), Pd{k) and Ptx{k) the respective 
probabilities that two pairs of (x, y) independently drawn from Zk and Zq with replacement 
are concordant, disconcordant and tie in x at the fc-th measurement time. The serial Pk 
index at the k-th measurement time is denoted as 


Pk(A:) 


Pc{k) + ^Ptx{k) 

Pc{k) + Pd{k) + PUk)' 


The PK(fc) value is interpreted in the same way as that in the Pk analysis. A value of 
one means that the indicator always correctly predicts the observed depth of anesthesia, 
a value of 0.5 means that the indicator predicts no better than 50/50 chance, and a Pk(^) 
value less than 0.5 means that the indicator predicts inversely. To be more precise, a zero 
value also means a perfect prediction but in a reversed direction. 


y%)} 


The estimation of the serial Pk and its standard error was obtained by the Jackknife 
method. Before the null hypothesis, the serial Pk value less than 0.5 was converted into 
one minus the serial Pk value. 
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Results from Clinical Database 


6.1 Visual Information in TF Plane 

Abundant information can be visually read from the tvPS of IHR. Various features were 
determined to reflect the physiologic dynamics of human body responding to various anes¬ 
thetic events (Fig. 6.1). 

It is demonstrated (Fig. 6.2,6.3) that the tvPS is more concentrated on deeper anesthetic 
levels and more scattered on lighter anesthetic levels. This then created the hypothesis 
on trying to understand the “rhythmic-to-non-rhythmic” phenomenon. This phenomenon 
is consistently seen in most study cases. During spontaneous breathing (Fig.6.3), the 
dominant curves can be seen in the first portion (usually on the left hand side) of the tvPS 
graph. Meanwhile, less dominant curves can be identihed in the second portion (usually 
the right hand side). The dominant curve represents the typical rhythmicity feature. 
It is also important to observe that the tvPS becomes more non-rhythmic prior to the 
appearance of first reaction. (Fig. 6.3) 

Painful surgical stimulations during anesthesia are referred to as noxious stimulation be¬ 
cause an unconscious patient under anesthesia cannot report pain sensation. Two of the 
most intense noxious stimulationthat occur in anesthesia and surgery are endotracheal in¬ 
tubation and skin incision. Clinical anesthesiologist uses the increase of heart rate and 
blood pressure as indicators for the level of noxious stimulation. In my study, we have 
found the visual information of tvPS that reveals the influence of noxious information 
where the low frequency (LF) power increases immediately after noxious stimulation, and 
wanes soon. The “TF surge phenomenon'’^ appeared in all study cases. Although the main 
goal of my study is on investigating NRR index, it can be seen that the potential value of 
LF power as an “noxious index” deserves to be investigated too. 
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Painful stimulus Time-varying Sharp change in time 

multi-components 


Figure 6.1: The time-varying spectrum of R-R peak interval (RRI) recorded during 
anesthesia showing rich dynamic features of human body in anesthesia. The RRI tracing 

is superimposed as red line 


These visual findings are quantified by NRR and HRV indices, such as tvHF, tvLF and 
tvLHR for subsequent analyses. 


6.2 Multiple Component Phenomenon 

Occasionally, multiple independent oscillatory components hidden in IHR can be seen 
from the tvPS. To further expand on this statement, the term independent is defined 
as the frequency of one component of which the modulation is unrelated to other com¬ 
ponents. Therefore, this is definitely not the harmonics (Fig.6.8, 6.9, 6.10, 6.11). The 
“multiple component phenomenon” are an interesting finding, which demonstrates the 
dynamic physiology of human body undergoing anesthesia. 

Apparently, it is nearly impossible to identify the appearance of multiple component from 
raw IHR waveform using naked eye. The tvPS is good at presenting “multiple compo¬ 
nent phenomenon”, but it is a difficult task for classical power spectrum. The “multiple 
component phenomenon” appeared in one thirds of the study cases. The moment that it 
appeared and disappeared seems random and unpredictable. Even so, this phenomenon 
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Figure 6.2: A representative data from a patient under anesthesia and controlled venti¬ 
lation. The 1200 seconds R-R interval (RRI) is superimposed on its time-varying power 
spectra (A) as a blue solid line. The corresponding electrocardiography derived respira¬ 
tion (EDR) is superimposed on its time-varying power spectra (B) as a green line. Panel 
C shows simultaneous recorded Bispectral index (BIS), effect-site sevoflurane concentra¬ 
tion (Sevo) and non-rhythmic to rhythmic ratio (NRR). Arrow 1: inadequate level of 
anesthesia was noted and the sevoflurane concentration was raised. Sevoflurane concen¬ 
tration reached top level after 300-second. Arrow 2: both BIS index and NRR reached 
minimum. There is a dominant straight line around 0.167 Hz in the panel A, whose 
fundamental frequency is the same as the rate of the ventilator, representing the en¬ 
trainment of RSA by the controlled ventilation. The time-varying power spectra of RRI 
showed rhythmic-to-non-rhythmic transition, corresponding to varying anesthetic depth, 
whereas the time-varying power spectra of EDR was almost unchanged under controlled 

ventilation. 
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Figure 6.3: A representative data obtained from the same patient under spontaneous 
breathing in emergence period. The 1200-second R-R interval (RRI) is superimposed on 
its time-varying power spectra (A) as a green solid line. The corresponding electrocar¬ 
diography derived respiration (EDR) is superimposed on its time-varying power spectra 
(B) as a blue solid line. Panel C shows simultaneous recorded Bispectral index (BIS), 
effect-site sevoflurane concentration (Sevo) and non-rhythmic to rhythmic ratio (NRR). 
Arrow I: endotracheal suction and extubation, Arrow 2: first reaction, Arrow 3: regain 
of consciousness. Under spontaneous breath, RRI and EDR are similar in high frequency 
region of spectrograms. The rhythmic component and its harmonics are time-varying in 
both time-varying power spectra of RRI and EDR. With the overall decrease of sevoflu¬ 
rane, both RRI and EDR exhibited a rhythmic to non-rhythmic transition. The RRI and 
EDR became non-rhythmic before the first reaction. (From ref. [69]) 























Chapter 6. Results 


37 



0 250 500 750 


0 250 500 750 

Time (s) 

-LF Power-LHR BIS 


Figure 6.4: A representative data from a patient under anesthesia receiving skin incision 
for her surgery. The 750-second R-R interval (RRI) is superimposed on its time-varying 
power spectra (A) as a green solid line. The corresponding electrocardiography derived 
respiration (EDR) is superimposed on its time-varying power spectra (B) as a blue solid 
line. Panel C shows simultaneous recorded Bispectral index (BIS), effect-site sevoflurane 
concentration (Sevo) and non-rhythmic to rhythmic ratio (NRR). 

is related to specific anesthetic events like only appearing when the patient is undergoing 
anesthesia. Although mechanism underlying this phenomenon is beyond the goal of my 
present study, it still presents an interesting research topic for the future. 

6.3 Ability to Predict Anesthetic Events 

In the serial Pk analysis for the first reaction (Fig. 6.14), the NRR is ahead of BIS 
{p < 0.05 for 20 seconds). The NRR was able to predict the first reaction 30 seconds 
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tvPS of IHR 



Figure 6.5: tvPS of IHR during emergence period from the last three consecutive cases 
(no.l) showing rhythmic-to-non-rhythmic phenomenon. BIS index (golden yellow tracing) 
and IHR (blue tracing) are superimposed. 


tvPS of IHR 



Figure 6.6: tvPS of IHR during emergence period from the last three consecutive cases 
(no.2) showing rhythmic-to-non-rhythmic phenomenon. BIS index (golden yellow tracing) 
and IHR (blue tracing) are superimposed. 
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tvPS of IHR 



Figure 6.7: tvPS of IHR during emergence period from the last three consecutive cases 
(no.3) showing rhythmic-to-non-rhythmic phenomenon. BIS index (golden yellow tracing) 
and IHR (blue tracing) are superimposed. 



Time (s) 

Figure 6.8: The time-varying power spectrum of instantaneous heart rate (IHR) 
recorded during anesthesia showing three independent oscillatory components, revealing 
rich dynamic features of human body in anesthesia. 
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Figure 6.9: The IHR (blue tracing) was recorded during spontaneous breathing. The 
second component appeared shortly (200-300 s). 


tvPS of IHR 



Figure 6.10: The IHR during spontaneous breathing showed “rhythmic-to-non- 
rhythmic” changes with BIS index (golden yellow tracing), accompanied with multiple 

components during “rhythmic” period. 
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tvPS of IHR 



Figure 6.11: Another IHR data (blue tracing) during spontaneous breathing showed 
“rhythmic-to-non-rhythmic” changes with BIS index (golden yellow tracing), accompanied 
with multiple components during “rhythmic” period. 

in advance (Pk > 0.90). At the instance of first reaction (0 seconds), both the NRR 
and BIS (both PK maxima > 0.95) were significantly better than other parameters. The 
time-varying HRV indices and HR (Pk < 0.83) provides significantly worse results than 
the NRR. 

In the serial Pk analysis for skin incision (Fig. 6.4), tvLF reaches maximum (Pk > 0.95) 
roughly at 30 seconds after skin incision. It reflects skin incision best, followed by tvLHR 
(Pk < 0.85). The tvLF is significantly better than BIS (Pk < 0.55) and NRR (Pk < 0.65). 

In the serial Pk analysis for LOG (Fig. 6.12), BIS reflects perfectly (Pk = 1) 50 seconds 
after LOG. LOG is also associated with a decrease in tvLF, tvHF and HR. Nowever, NRR 
does not reflects LOG well. In the serial Pk analysis for ROG (Fig. 6.15), BIS is the best 
index, and surpasses NRR, HR and the time-varying HRV indices significantly. 

Due to the inadequate data quality, 1 patient has been eliminated from our data in regards 
to the serial Pk analysis of the initial reaction. For our study, when we analyzed the serial 
Pk of skin incision, LOG and ROG, we also eliminated 3, 5, and 4 patients respectively. 


6.4 Correlation with Sevoflurane Concentration 

The correlation with the estimated effect-site sevoflurane concentration is evaluated for 
the spontaneous breathing (SB) period (table 6.2). The data number of SB is 5810 (23240 
s) and the average duration for each patient is 750 ± 322 s. The overall indices ranking of 
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Figure 6.12: Tracings and their faded error bars about the prediction probability (Pk) 
values of the parameters and their standard errors in serial Pk analysis for loss of con¬ 
sciousness (LOG). The baseline is one minute before LOG. NRR = non-rhythmic to 
rhythmic ratio; HF = high frequency power; LF = low frequency power; LHR = LF to 
HF ratio; BIS = Bispectral Index; HR = heart rate. (From ref. [69]) 


the Pk analysis and Speareman rank correlation (R) are consistent. While BIS correlates 
with effect-site sevoffurane concentration best (SB: Pk = 0.839, R = —0.836, p < 0.0001), 
the NRR index is second (Pk = 0.736, R = —0.619, p < 0.0001). The tvLF is the best 
HRV indices ( Pk = 0.669, r = —0.416, p < 0.0001). Compared to tvLF, the NRR index 
is significantly better ( p < 0.01). BIS surpasses NRR significantly in SB (p < 0.001). The 
NRR index is significantly better than HR (p < 0.0001). 


6.5 Quantitative Results of Noxious Stimulation 

When observing the serial Pk analysis for skin incision (Fig. 6.4), tvLF reaches the 
maximum amount (Pk > 0.95) 30 seconds after skin incision. This shows that the tvLF is 
the best index of the skin incision, followed by the heart rate (Pk < 0.85). Furthermore, 
tvLF is significantly better than BIS (Pk < 0.55) and the NRR index (Pk < 0.65). 
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Figure 6.13: Tracings and their faded error bars about the prediction probability (Pk) 
values of the parameters and their standard errors in serial Pk analysis for skin incision. 
The baseline is 5 seconds before skin incision. NRR = non-rhythmic to rhythmic ratio; 
HF = high frequency power; LF = low frequency power; LHR = LF to HF ratio; BIS = 
Bispectral Index; HR = heart rate. *P < 0.05, LF versus BIS and LF versus NRR. (From 

ref. [69]) 


In the serial Pk analysis for endotracheal intubation (Fig. 6.16), I compared different 
indices calculated from multitaper reassignment spectrogram (HF, LF, LHR) with the in¬ 
dices calculated from STFT spectrogram (HFstft j LFstft , LHRstft )• This analysis 
signifies that the LF power and LHR provide the best results for the noxious stimulation 
of endotracheal intubation. Meanwhile, the heart rate is represented the second best and 
the BIS index is represented the worst. The indices that are calculated from multita¬ 
per reassignment spectrogram performs better than the corresponding indicators that are 
calculated from the STFT spectrogram. 


6.6 Summary 


Several features on TF plane can be easily seen conveniently via the tvPS, including 
the “rhythmic-to-non-rhythmic phenomenon”, “multiple component phenomenon”, and 
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Figure 6.14: Tracings and their faded error bars about the prediction probability (Pk) 
values of the parameters and their standard errors in serial Pk analysis for first motor 
movement reaction. The baseline is 3 minutes before first motor movement reaction. 
NRR = non-rhythmic to rhythmic ratio; HF = high frequency power; LF = low frequency 
power; LHR = LF to HF ratio; BIS = Bispectral Index; HR = heart rate. * P < 0.05, 

NRR versus BIS. (From ref. [69]) 


“LF surge phenomenon”. The “rhythmic-to-non-rhythmic phenomenon” and “LF surge 
phenomenon” are apparent in the majority of my study cases. This universal consistency 
foreshadows the subsequent positive statistical results. 

The quantitative analysis proves that the NRR index indicates the first motor movement 
reaction better than other clinical indices [69]. BIS index is the best indiator for LOG 
and ROC. NRR index correlates with the concentration of sevoflurane. Furthermore, LF 
power indicates noxious stimulation best, and better than the heart rate. It is also shown 
the statistical benefit of the multitaper reassignment spectrogram over the classical STFT 
spectrogram. 
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Figure 6.15: Tracings and their faded error bars about the prediction probability (Pk) 
values of the parameters and their standard errors in serial Pk analysis for return of 
consciousness (ROC). The baseline is 3 minutes before ROC. NRR = non-rhythmic to 
rhythmic ratio; HF = high frequency power; LF = low frequency power; LHR = LF to 
HF ratio; BIS = Bispectral Index; HR = heart rate. * P < 0.05, BIS versus NRR. (From 

ref. [69]) 
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Figure 6.16: Tracings and their faded error bars about the prediction probability (Pk) 
values of the parameters and their standard errors in serial Pk analysis for endotra¬ 
cheal intubation. The baseline is 5 seconds before endotracheal intubation. HF = high 
frequency power calculated by multitaper reassigned spectrogram; LF = low frequency 
power calculated by multitaper reassigned spectrogram; LHR = LF to HF ratio calcu¬ 
lated by multitaper reassigned spectrogram; HF_stft = high frequency power calculated 
by STFT spectrogram; LF_stft = low frequency power calculated by STFT spectrogram; 

BIS = Bispectral Index; HR = heart rate. 
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Table 6.1; Representative values of the Indices during important anesthetic events. 
Values are presented as median values (lower and upper quartiles). HF, LF, and LHR are 
expressed in common logarithm (10 base) of millisecond square. Heart rate is expressed 
in beat per minute. BIS = Bispectral Index; NRR = non-rhythmic to rhythmic ratio; HF 
= high frequency power; LF = low frequency power; LHR = LF to HF ratio; LOG = loss 
of consciousness; ROC = return of consciousness; T0= 60 s before LOG; Tl= 60 s after 
LOG; T2= 5 s before incision; T3= 30s after incision; T4= 180 s before first reaction; 
T5= 5 s after first reaction; T6= 180 s before ROC; T7= 5 s after ROC. (From ref.[69]) 
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0.79 

0.68 

0.91 


(1.50, 

(1.01, 

(0.46, 

(1.25, 

(-0.37, 

(0.75, 

(0.01, 

(0.37, 


2.27) 

1.60) 

0.94) 

2.25) 

0.56) 

1.31) 

1.27) 

1.26) 

LHR 

-0.42 

-0.03 

-0.50 

0.24 

-1.33 

-0.25 

-0.44 

-0.34 


(-0.84, 

(-0.29, 

(0.46, 

(0.01, 

(-1.67, 

(-0.73, 

(-0.79, 

(-0.66, 


-0.27) 

0.19) 

0.84) 

0.50) 

-0.66) 

0.06) 

0.10) 

-0.12) 

HR 

74 (65, 

64 (58, 

65 (58, 

70 (60, 

71 (62, 

79 (69, 

76 (69, 

84 (75, 


91) 

72) 

72) 

76) 

80) 

85) 

83) 

91) 
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Table 6.2: Prediction Probability and Spearman Rank Correlation for Indices vs. Es¬ 
timated Sevoflurane Effect-site Concentration in Controlled Ventilation and Spontaneous 
Breathing. * p < 0.05 between non-rhythmic to rhythmic ratio (NRR) and other indices. 
Indices with prediction probability (Pk) value smaller than 0.5 was corrected by 1 — Pk 
before comparison. BIS = Bispectral Index; tvHF = high frequency power; tvLE = low 
frequency power; tvLHR = LF to HE ratio; HR = heart rate. (From ref. [69]) 


Controlled Ventilation 


Pk (SE) 

Spearman rank correla¬ 
tion (95% Cl) 

BIS 

0.716 (0.020) 

-0.575 (-0.661, -0.476) 

NRR 

0.670 (0.025) 

-0.467 (-0.583, -0.337) 

tvHF 

0.479 (0.027)* 

0.073 (-0.069, 0.211)* 

tvLF 

0.582 (0.025) 

-0.233 (-0.359, -0.101)* 

tvLHR 

0.581 (0.025) 

-0.233 (-0.364, -0.096)* 

HR 

0.423 (0.024)* 

0.152 (0.034, 0.269)* 


Spontaneous 

Breathing 


Pk (SE) 

Spearman rank correla¬ 
tion (95% Cl) 

BIS index 

0.839 (0.014)* 

-0.836 (-0.881, -0.778)* 

NRR 

0.736 (0.019) 

-0.619 (-0.688, -0.536) 

tvHF 

0.489 (0.023)* 

0.001 (-0.117, 0.113)* 

tvLF 

0.669 (0.022)* 

-0.416 (-0.506, -0.317)* 

tvLHR 

0.666 (0.022)* 

-0.420 (-0.514, -0.319)* 

HR 

0.532 (0.022)* 

0.027 (-0.075, 0.133)* 



Chapter 7 


Implications to Anesthesiology 
and Physiology 

7.1 Main Findings 

There are two important clinical results in regards to the NRR index. First, the NRR 
index is able to predict first reaction during emergency period. Second, the NRR index 
correlates with the sevoflurance concentration during the patient’s spontaneous breathing. 
Even so, the NRR index is not the best predictor for LOG and ROC. These findings suggest 
that there is a distinct physiological interpretation of the anesthetic depth level when 
compared to the hypnosis measured by surface EEG. Furthermore, tvLF correlates well 
with skin incision and endotracheal intubation. This momentary information is currently 
unavailable to clinical anesthesiologists as both the IHR and tvPS are not implemented in 
current monitor. 


7.2 Genuineness of NRR index 

The NRR index generally has a sharp increase nearly before the onset of the patient’s 
first reaction. This then makes it highly unlikely that the motion artifacts interfere the 
NRR index. Furthermore, the observed NRR index performance in not confounded by 
the ventilation mode change. In particular, aside from the LOG, there are no transitions 
that occur from the mechanical ventilation to the spontaneous breathing in the data we 
have obtained. Also, there is no positive pressure ventilation for all the first reaction 
data. Furthermore, past literature on the association between ” rhythmic-to-non-rhythmic” 
transition of the oscillatory pattern in RRI index and the different level of anesthesia has 
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been mentioned but no quantitative works have been done. By using the classical power 
spectral analysis, Kato et al. has described that the high frequency component of RRI 
disperses when the patient awakens and concentrates when the patient is under anesthetics 
and controlled ventilation[2], Based on the above facts, I am more than confident to state 
that the NRR is convincing in its distinct and novel information that it is capable of 
providing us. 


7.3 Existing Findings 

The results we have obtained from the BIS (Fig.6.12, Table 6.2) are in agreement with 
previous studies showing that BIS is capable of monitoring awake status vs. anesthesia, 
the first reaction and the decrease of anesthetic gas concentration[5, 8, 10]. Our BIS results 
(Fig.6.13, Fig.6.16) also agree with past studies that provide evidience that BIS or other 
EEG-derived indices cannot indicate the response to noxious stimulation [6, 7, 9]. 

The results from our analysis (Table 6.1) are also consistent with previous reported associ¬ 
ation between anesthetic depth and traditional HRV parameters: HRV is decreased during 
general anesthesia and increased during recovery[2, 46, 70, 71]. Lastly, our results are also 
consistent with the reported finding that combining classical HRV indices with BIS does 
not outperform BIS in discriminating awake from asleep during anesthetic induction[72]. 

It has been reported that classical HRV parameters do not outperform heart rate when 
predicting noxious stimulation[4, 73]. Contrarily, our results demonstrate that tvLF pre¬ 
dicts skin incision better than HR. It seems that the noxious stimuli of the skin incision 
elicited a momentary increase of sympathetic activity, which leads to this finding. Al¬ 
though the study is not designed for noxious stimulation, the finding that tvLE correlates 
with skin incision shows the technical advantage of our approach and its potential as an 
index for noxious stimulation. 


7.4 Physiologic Interpretation of NRR 

The association between NRR and anesthetic depth can be partially explained by the 
cardiopulmonary coupling[14] and the respiratory physiology. When the coupling effect 
is higher, the respiration pattern is reflected in the RRI index. The respiratory mecha¬ 
nism also contributes to the differential influence that occur from anesthetics. The neural 
respiratory control comprises of two systems, the involuntary automatic control system 
and the voluntary control system. The involuntary automatic control system is mainly 
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controlled by the respiratory center in the pontomedullary area. The voluntary control sys¬ 
tem is mainly controlled by the forebrain [25]. These two systems are also distinct in their 
neural anatomy and their functions. Respiratory efferent signals from these two systems 
compete with each other and are integrated at the spinal level to control the respiratory 
motor neuron. However, during spontaneous respiration, studies have also shown that 
the involuntary automatic respiratory pattern generated in the pontomedullary center [25], 
specifically the preBotzinger complex[26-28], is rhythmic. 

On the other hand, the breathing pattern of the voluntary respiratory motor control is 
non-rhythmic and involves the cortical processing and thalamic integration which responds 
to the peripheral inputs and the descending inputs [25, 30-32]. The thalamus also ac¬ 
tively participates in different motor functions, including respiration[31, 32], speech[33] 
and cough[34] to name a few[74]. Our NRR index findings suggest that the non-rhythmic 
respiratory activity involves the nearly entire cerebral regions is more susceptible to sevoflu- 
rane than the rhythmic respiration generated in the medulla. This suggestion is supported 
by literature. Guedel made the classical observation on respiratory pattern[l]. Bimer et 
al denoted that the respiratory irregularities accompany electroencephalographic arousal 
reaction in spontaneous breathing[75]. Also, Studies have also shown that the respiratory 
activity of preBotzinger complex is less depressed by sevoflurane[36, 37], and the effect of 
volatile anesthetics is less prominent to the brainstem compared to the cortical regions 
that are more abundant in synaptic transmissions. [37]. 

The ability for the NRR index to predict first reaction suggests a connection between the 
“rhythmic-to-non-rhythmic” phenomenon and motor reactions. Antognini et al. demon¬ 
strated that immobility to surgical stimulation by volatile anesthetics in goats is mainly 
modulated by the spinal cord, which also indirectly affects the thalamic response to nox¬ 
ious stimulation[12, 13]. Another study by Velly et al. also presented significant results[7]. 
He recorded the human subcortical electrophysiological activity and demonstrated that 
subcortical activity (possible the thalamic activicy) predicts suppression of movement to 
noxious stimuli, but not changes of consciousness, whereas cortical EEG predicts loss 
of consciousness, but not motor suppression. Based on the above, NRR might reflect 
the subcortical activity since it better reveals information about motor reaction, but not 
consciousness, compared with BIS. Since immobility to noxious stimulation is similar to, 
although not the same as the first reaction under surgical wound pain in the present study, 
this relationship suggests the possible role of the thalamus in “rhythmic-to-non-rhythmic” 
phenomenon. Besides the NRR index, BIS also responded to first reaction in our serial 
Pk analysis, which is in agreement with other studies[10]. Therefore, it is possible that a 
cortical mechanism is also related to the first reaction. 
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From these results, it is possible to see that sevoflurane may affect non-rhythmic respiration 
level of the spinal cord, subcortical supraspinal regions, or the cortex in the human brain. 
Although we are currently unable to pinpoint the exact anatomic location, it is likely that 
the subcortial supraspinal area, possibly the thalamus, is the most plausible region. The 
above statement is sound as the NRR index generates different results in the serial Pk 
analysis and is also closely connected to the human body’s respiration responses. Thus, we 
hypothesize that the “rhythmic-to-non-rhythmic” phenomenon in IHR exhibits the central 
respiratory activity via either central or peripheral mechanisms [14]. The IHR also reflects 
the final integration of this ” rhythmic-to-non-rhythmic” respiratory activity. Although 
more evidence is necessary to further clarify these hypotheses, we propose to use the NRR 
quantification methodology as a potential tool to evaluate the depth of anesthesia that 
differs from EEG-based monitoring. 


7.5 Clinical Application of the NRR Index 

The first strength of this study is a new quantitative approach, referred to as the NRR 
index, to analyze “rhythmic-to-non-rhythmic” phenomenon. The NRR index has the po¬ 
tential to reflect different levels of anesthesia. Literature partially supports that there 
is an underlying physiological mechanism. Second, the signal processing technique, mul¬ 
titaper Synchrosqueezing transform, has been well studied in the literature, so we have 
adequate theoretical support to resolve the potential limitations of the PS. The momentary 
dynamics in IHR can be captured by the tvPS, which leads to the NRR index. 


7.6 Potential Index in Noxious Stimulation 

Although we did not design the study for noxious stimulation, from the data analysis 
result, we found that the indices we apply are potential in detecting pain reaction under 
anesthesia. In particular, the tvLF index correlates with the noxious stimulation better 
than heart rate as an additional finding. More study will need to be conducted on the 
tvLE to provide sound results. It can also be helpful to compare this data with other 
well known indices like the Surgical Stress Index® (SSI, GE healthcare) and Analgesia 
Nociception Index® (ANI, Metrodoloris), which have been used on monitoring surgical 
stress [11]. SSI is derived from the pulse wave amplitude and pulse beat interval of the 
photoplethysmography whereas ANI indicates the parasympathetic tone derived from EGG 
signal. It is possible that various indices can provide complementary information regarding 
noxious stimulation that is gathered from different modalities and concepts. 
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NRR index measures the rhythmic-to-non-rhythmic transition and quantifies a new kind of 
information different from the above two instruments in both physiologic and mathematical 
senses. Further investigations may be necessary to understand the clinical value and 
application of tvPS. 


7.7 Summary 

In conclusion, we are able to extract hidden information regarding the anesthetic levels 
from the routine ECG monitor. The quantitative results of the NRR index supports our 
initial hypothesis that the “rhythmic-to-non-rhythmic” transition correlates with motor 
reactions during emergence period earlier than BIS index, and correlates with sevoflu- 
rance concentration. In addition, the potential of the time-varying HRV indices provides 
observations of the relationship between tvLF index and the noxious stimulation. The 
notion of these dynamics is rooted in the recently developed signal processing technique. 
Without tvPS, the above quantification cannot be easily achieved. Overall, my clinical 
study suggests that the ECG signal contains complementary information to the EEC-based 
depth-of-anesthesia index. 



Chapter 8 


Real-time Processing Using the 
Blending Operator 

8.1 Real-time Processing 

From the previous chapters, my study has concluded that the IHR signal during anesthesia 
contains a wealth of information when viewing with tvPS. It is possible to derive the 
IHR signal from ECG, the waveform of pulse oximetry, and also the invasive intra-arterial 
pressure waveform. Although we mainly use ECG signal in the present project, all possible 
sources of IHR are mandatory or frequently used functions in modern standard anesthesia 
monitors. A practical question is how do we apply the NRR index proposed in the present 
study to individual patients. 

The processing steps of NRR index starts from the identification of each heart beat from 
ECG signal. Next, an interpolation is used to convert the irregularly-sampled heart beats 
data into regularly-sampled IHR data. Next, the calculation of the tvPS is based on the 
multitaper Synchrosqueezing tranform. And lastly the NRR index is computed based on 
the tvPS. Thus, the positive statistical results in previous chapter, are all obtained from 
off-line calculations. 

In anesthesia practices, the main purpose of anesthesia monitor is to provide real time 
monitoring of the patient’s responses in order to provide timely responses. This purpose 
then turns ordinary measurements like heart rate, blood oxygen saturation, and blood 
pressure into vital information during anesthesia. The real-time processing of NRR in¬ 
dex is necessary if we would like to “steer” the anesthetic depth of patients based on 
their instantaneous information from IHR. The potential benefits include adequate dose 
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of anesthetics, possibility to divert critical situation, reduce stress response, and improve 
long-term welfare. 

In this chapter, my main focus is the real-time processing from the irregular heart beat 
sample to the tvPS. The goal is to interpolate the irregular data samples in real time. 
Meanwhile, the polynomials property of the interpolation also serves as the vanishing- 
moment and minimum-supported wavelets. 


8.2 Technique Obstacles for Real-time Processing 

Due to the irregular interval of each heartbeat, the heartbeat data is irregularly sampled 
and this irregularity is unavoidable. Almost all continuous information are digitized in 
equally-spaced samples in time or in space before further analysis for convenience, and for 
efficiency. However, some situations provide non-uniform sampling and could be either 
unavoidable or intended[76, 77]. In order to derive the continuous instantaneous heart 
rate from the distinct heartbeats (Fig.2.1) is an example of the rich IHR information that 
entirely relies on the irregularity of the heartbeats. Hence, the initial challenge is the 
real-time processing of the non-uniform data samples. 

When analyzing data the scope of continuous wavelet transform (CWT), the theoretical 
background of Synchrosqueezing transform presents the optimal performance when us¬ 
ing compactly supported windows in frequency domain [44]. That means we need whole 
data converted in frequency domain for subsequent calculation. This is inefficient and 
also highly unlikely for real-time data processing. Moreover, for real-time processing, an 
wavelet transform (or STFT) better has a short, and compact support window in time. 
This type of window may hamper the benefit of wavelet transform; an adequate vanishing 
moment can eliminate trend component in terms of corresponding polynomial order to 
preserve the oscillatory components we desire more efficiently. 

A possible workaround is by cutting down a segment of data, windowing it, and then 
choosing a suboptimal wavelet function empirically. Next, we would observe the effect of 
suboptimal window and then adjust the parameters empirically until we finally accept the 
suboptimal effect in real-time display. 

There is a better solution to consider before we proceed with the above workaround. The 
solution is a special spline interpolation, which has the interpolating ability to handle the 
unequally-spaced heart beat sample, and the polynomial property to match the vanishing 
moment property of wavelet transform. The resulting time-scale analysis is compactly 
supported in time and provides better performance due to the vanishing moment property. 
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Based on this scheme, the real-time tvPS can help remove the above mentioned technical 
obstacles. 


8.3 B-spline 


We start from a brief introduction of the spline function. For digital samples g{tn) are 
sampled at time points 

t : a = to < h < t 2 < ■ ■ ■, (8.1) 

where {t^} may be irregular, or non-uniform signal, which means that we allow t^+i —t^^ 
tj+i — tj for k ^ j. Let m > 3 be any desired integer and 11^-1 denote the space of all 
polynomials of degree less than m. 

The spline space is then defined as: 


s ; S-m+l = S-m+2 = ■ ■ ■ = a = So < Si < S2 < ■ ■ ■, 


( 8 . 2 ) 


and the notation Ss,m ■= <S's,m[a, co) denotes the spline space of order m with the knot 
sequence s. 

To formulate a locally supported basis of Ss,mi we now introduce the truncated powers 

:= (max{0,x})™-^ 


The truncated power functions (sfc — t)+~^, /c = 0,1,..., are in Ss,m- Since {x^ — 

A: > 1, have global support, we apply the m-th order divided differences to form locally 
supported functions. The divided differences are 


[u,...,u]/ := 


ll 


if there are I -I- 1 entries in [u,... ,u], and 


[^^0, 


[ui, ...,Un]f- [no, . . . ,Un-l]/ 
Un - Uo 


if uo < ui <...< Un and Un > uq. 


Now we consider the knot sequence s of the spline space Ss,m- The normalized B-splines\7S\ 
are the divided difference of the truncated power (s^ — namely. 


d^m,kif ) (•Sm+fc 'Sfc) [ Sfc, . . . , ] (• t)_|_ , (S.3) 


for k = —m -|- 1,..., 0,1, 2,.... 
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We would apply the recursive formula to implement Nm,k [52, page 143 (6.6.12a)]. The 
B-splines are a locally supported. 

However, if we choose the irregularly-spaced sequence t of time positions in (8.1) as the 
knot sequences s in (8.2) by attaching t-m+i = ... = to = a to t as in (8.2), it can be 
difficult to compute, or the computation cost will be higher for the non-uniform {tj} and 
large values of n. 


8.4 Quasi-interpolation 


On the other hand, de Boor and Fix [79] proposed a spline representer f{t), which ap¬ 
proximates the data instead of providing exact interpolating the target data function g{t) 
at t = tj, j = 0,1,... ,n (that is, if f{tj) / g{tj) is allowed). This quasi-interpolation has 
a “polynomial preservation” property; g{t) = p{t) £ n^-i. 

This polynomial preservation is important to improve the performance of the continuous 
wavelet transform (CWT), and hence the Synchrosqueezing transform (SST), to reveal 
the oscillatory components. Specifically, if the analysis wavelet is compactly supported by 
the spline of order m, then the m-th order vanishing moment annihilates the m-th order 
Taylor polynomial expansion of g{t). This then facilitates the trend removal of the the 
signal g{t) as a polynomial. However, the quasi-interpolation scheme de Boor proposed 
[78, 79] requires derivative data values of g{t). Others [80, 81] proposed to replace the 
derivatives of g{t) by divided differences of {g{ti)}. However, these quasi-interpolations[82] 
do not address the real-time issue at stake. In addition, quasi-interpolation introduced the 
error g{ti) - f{ti) for i = 0,1,..., n. 

Chui et al. proposed a real-time quasi-interpolation [83]. Chui and his college also intro¬ 
duced a “local interpolation” [84] to correct the error of “quasi-interpolation”. 

The real-time quasi-interpolant [83] is introduced as follow; for each A: = 0,1,..., consider 
the Vandermonde determinant 


• 1 tk+m+l) • det 


1 

tk 


1 

tk+l 


1 

tk+m—1 


j.m—1 j.m—1 

^k ''fc-l-1 


t 


m—1 
k+m—1 J 


and the determinant D{tk, ■ ■ ■, tk+j-i,Cj,tk+j+i,... ,tk + rn — l) obtained by replacing the 
(j -|- l)-st column in the definition of D{tk, ■ ■ ■, tk+m+i) by the column vector 
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where = 1 and 


Cij,rn) 


<7 {tj+li ■ • • tj+m—l) 



for i = 1,... ,m — 1 with cj*(ri,..., rm-i) being the classical symmetric functions defined 
by , Tm-i) = 1 and for i = 1,..., m — 1, 




Using the determinants introduced above, we apply the below spline coefficients 

_ ‘ ' 1 ^k+j+l-) • • * ; l) 

D{tk, ■ ■ ■ ,tk+m-l) 

to formulate the compactly supported spline function 

2m—2 

(*):= E 0-k,l—m+lRlt,m,k+l—m+l {t) i 

l=m—l 

with suppMf^m,fc = [tk-m+i-,tk+m\- These basis functions provide a real-time implementa¬ 
tion of the quasi-interpolation operator 

(Qmff)(t) = y^^g(tk)Mt,m,k(t). (8.4) 

k 


The quasi-interpolation operator Qm (8.4) possesses the polynomial preservation property. 


For any m > 1, 


(QmP)(t) =p(t), 


for all t > a and for all p G Um-i, provided that the summation (8.4) account for all non¬ 
negative integers k = 0,1,.... Furthermore, in view of the support of Mt^m,k-, h follows 
that 

v-\-ra—2 

P{tk)Mt,m,k{t) = pit), tG[tu,t.u] 

k=v—m-\-l 

for all p G Ilm-i. 


This local polynomial preservation property allows the CWT, and hence the SST, to 
annihilate the (m — l)-th degree Taylor polynomial approximation of the signal at tj, 
where u < j < v. 
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8.5 Local Interpolation Operator 


The local interpolation operator, denoted by Rm, satisfies the interpolation property. To 
define Rm, we will insert knots to t by considering a new knot sequence sDt (8.2). For 
simplicity, we will introduce even order variable m : 

For even m > 4, we set 

•Smfc/2 k 0,1, 2,.... 

That is, we insert (m/2 — 1) knots in between two consecutive knots in t. For convenience, 
we will choose the new knots to be equally spaced between every pair of two consecutive 
knots. 

Fix even m > 4. Let be the m-th order B-spline with knot sequence s. Then the 

completely local spline basis function can be defined by 

T u\ — 

s,m,m{j—l) /2IL / 

Since tj = s^j /2 is the “centered” knot and 

SyippNs,m,m{j—l)/2 [1)/2) ^m(j + l)/2 ] — 


we have 

/ 

k-‘s,m,j (tj ) — 1 

< (8.5) 

SUppL/g^mJ — [tj — l:tjj^l]. 

It is clear that Ls,m,j(tk) = Sj-k, where the Kronecker delta notation is used (8.5). 

The above preparation provides a real-time implementation of the local interpolation op¬ 
erator, which satisfies the interpolation property (8.5). 

Fix m > 3. For a given function g G C'(M), the local interpolation operator R^ is defined 
as 

(^mg){t) := y^^g(tk)Ls,7n,k{t). ( 8 . 6 ) 

k 


8.6 Blending Operator 


We are now ready to apply (8.4) to obtain the blending operator, denoted by R^ © Qm- 
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Fix m > 3 and g G C'(M). The blending operator is defined as := Pm © Qm, where 
Pm © Qm •— Qm © Pm(l Qm) — Qm © Pm PmQm; 
and I is the identity operator. In particular, we have 

g{tk)Mt^rn,k{t) © Z - Z 

k k j 

We remark that in the definition of Pm, the two operators Pm and Qm are not commutative. 
Let us summarize the two key properties of the blending operator in the following theorem. 

1. The blending operator Pm possesses both the polynomial preservation property of 
Qm and the interpolatory property of Pm- 

2. The error of spline interpolation by the blending operator is both small and bounded. 


In conclusion, the blending operator, as a local spline interpolation operator, achieves the 
optimal interpolation error rate when compared with the traditional spline interpolation 
operator. In addition, the error depends only on the local data profile, which allows 
real-time implementation with optimal error rate. 


8.7 Real-time tvPS 

We provide a real-time computational algorithm to compute f{t) for the upcoming data 
samples g{to), g{ti), ■ ■ ■ etc. For convenience, we only consider even order m. The formu¬ 
lation for the odd order is similar but slightly more complicated. 

The blending operator leads to the real-time spline interpolation. For real-time issue, the 
wavelets have to be minimum supported in time. We can design the wavelets with the 
same number of vanishing moments as the spline. Meanwhile the wavelets have minimum 
support and maximum order of vanishing moments. 

Let m,n > I be arbitrary integers, and x an arbitrary knot sequence. The spline basis 
functions is then as follows: 




x^m+n^k 


(x), k G Z, 


( 8 . 8 ) 
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Algorithm 1 Real-time implementation of the blending operator 

First, pre-compute the B-spline values ='■ nij. Then for each k, since Mt,m,k £ 

Ss,m-, there exists a finite sequence {hk^i} in the formulation of 

I 


Also pre-compute 

Now, while the data sequence {^(tfc)} is acquired, compute 

9i = '^h,ig{tk) 

k 


and simultaneously compute 

* 9{ti) -Ejdijgitj) 

9i = -^-> 

^m(Z—1)/2) ^ 

and then up-sample {g^} by m{l — l)/2; that is, set 

9t(k-i)/2 = 9l and gf = 0. 

Then, we have an on-line computational scheme for the quasi-interpolation spline inter¬ 
polation: 

n 

f{t)= X] igi-9t)d^s,m,l{t) 

l=—m,+l 

for increasing number of samples from g{tn) to g{tn+i), ■ ■ •; and this can be implemented 
for real-time D/A conversion. 


to be called VM wavelets on x, where Nx,m+n,k{x) is defined in (8.3), which satisfies the 
moment conditions 


x''ipx,m-,n,k{x)dx = 0, / = 0,..., n - 1 


X 'tpx^m-,n,k(^x')dx ^ 0. 


(8.9) 


That means the arbitrary chosen n in the interpolator leads to the vanishing moment order 
of the wavelets if we use this spline basis as the wavelet function. The wavelet transform 
can successfully eliminate the n order polynomial as a trend component. 

Let m, n > 1 be arbitrary integers. Then the derivative of the wavelet '0£c,m;n,fc on an 
arbitrary knot sequence x, is given by 


'4’x,m-,n,ki^) ~ '4’x,m—l;n+l,k{x)- 


(8.10) 
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This property permits us to conveniently calculate the derivative of the wavelet, which 
is required to calculate the frequency reassignment in Synchrosqueezing transform, (see 
equation 4.3) 

Then, we can proceed to the Synchrosqueezed CWT to obtain the real-time time-frequency 
analysis as real-time tvPS. Using this tvPS, we can apply the research in previous chap¬ 
ters to quantify the “rhythmic-to-non-rhythmic” phenomenon in real-time. Using the 
clinical data to “emulate” the real-time calculation of NRR index, we obtained a compa¬ 
rable performance of the correlation with sevoflurane concentration during spontaneous 
respiration: The Pk value of real-time NRR index is 0.711 ± 0.021 (Pk of off-line NRR 
index:0.732±0.018)(6.2). 



Chapter 9 


Conclusion 


9.1 Research Findings 

Thanks to the advancement of science in anesthesiology and neural physiology, I can pro¬ 
pose a theory and modeling to investigate the “rhythmic-to-non-rhythmic” phenomenon. 
Furthermore, because of the advancement of signal processing technique, I am able to use 
multitaper Synchrosqueezing technique to compute the NRR index. Anesthesia creates 
dynamic changes to the human body, and due to the time-varying ability and good TF 
resolution, the NRR index is a particularly suitable tool reflecting the constant changes 
during anesthesia[43, 61]. 

The numeric result from clinical database is positive [69]. NRR index performs better and 
faster than BIS index when detecting motor movement reaction. This also supports the 
fact that different brain regions mediating various functions, with differential susceptibility 
to anesthetics may require multi-modal assessment of these anesthetic effects. NRR index 
is therefore able to provide additional value. In addition, LF index provides good results 
as an indicator for noxious stimulation, for both skin incision and endotracheal intubation. 
ECG waveform is standard information in current anesthesia, while the IHR is not. My 
study shows that these new indices derived from ECG waveform can provide additional 
information that could be beneficial to patients undergoing anesthesia in the future. 


9.2 Accomplishments 

1. The adaptive harmonic model is a method of modeling the rhythmicity of the “rhythmic- 
to-non-rhythmic” phenomenon 
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Figure 9.1: The diagram represents the overall framework from the human brain to 

signal analysis. 

2. A time-varying power spectrum that is based on multitapering Synchrosqueezed 
spectrogram. 

3. Methodology of the NRR index 

4. A methodology of serial Pk analysis was developed to evaluate the dynamic perfor¬ 
mance of anesthetic depth index as time-varying Pk value. 

5. Positive clinical value of the NRR index was obtained from a clinical database 

6. Improved performance of classical HRV parameters was obtained using tvPS 

7. A real-time scheme that provides real-time tvPS and real-time NRR index. 

9.3 Future Directions 

The results of my project indicates several future research possibilities. 

In clinical anesthesiology, more studies should be conducted to answer questions regard¬ 
ing the value of NRR index in different kinds of anesthetics (i.e. propofol, ketamine, 
dexmedetomidine or xenon), in different types of patient (i.e. children, lying-in women. 
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elderly patient, or patient with severe illness). Since the NRR index may indicate the 
activity of subcortical areas, this poses questions on if the application of the NRR index 
can lead to less stress response and better long-term outcomes, and if it is possible, how 
can we obtain those results. On the other hand, the IHR is derived from the ECG signal, 
I wonder whether it is possible to calculate the NRR index from photoplethsmography or 
intra-arterial pressure waveform. It also needs more clinical study to better understand 
the clinical role of NRR index and the tvLF index. 

It is possible that a thorough understanding of the basic mechanism can lead to better 
clinical benefit. The origin of “rhythmic-to-non-rhythmic” phenomenon should be in the 
brain. However, the neural mechanism has not been addressed according to my best knowl¬ 
edge. For example, where exactly is the anatomical location that creates this phenomenon? 
What could cause the multiple component phenomenon as previously mentioned? Lastly, 
how do the neurons coordinate with one another during this phenomenon? It is foreseeable 
that the answers of these questions can directly contribute to the brain science. Thus, the 
neurophysiological study based on the animal model is a potential future research. 

The signal processing technique directly contributes to the performance of these clinical 
indices. The current result of tvPS is helpful in revealing various features in the TF plane, 
but there are still a lot of room for improvements. The time-frequency analysis is based 
mostly on the Fourier transform. I am interested in whether it is possible to discard the 
idea of frequency in order to have a better understanding of the signal’s features. The real¬ 
time calculation and display of the clinical indices provides benefit to individual patients. 
Technique developments in real-time processing is imperative. 


9.4 Summary 

My study focuses on the qualifications of the “rhythmic-to-non-rhythmic” phenomenon. 
I have verified that my proposed model and quantitative index through using clinical 
database. The positive results support my proposal on using the adaptive harmonic model 
and time-varying power spectrum based on the multitaper Synchrosqueezing transform. 
The NRR index also presents unique values in clinical anesthesia. The achievements in 
this thesis fill a gap in the fields of medicine and provide a new perspective on many of 
the issues that could be understood in the future. 
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